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This report covers the period April 1 1991 to 
December 31 1991. Effort has continued through this period 
to refine and expand the SIRIUS/ABACUS program package for 
CASSCF and RASSCF second derivatives. One of the other 
authors of the package, Trygve Helgaker, and myself have 
devised a new approach to computing the Gaussian integral 
derivatives that require much of the time in gradient and 
Hessian calculations. This work has been written up and 
accepted for publication (copies attached) . 

Several different studies have been undertaken in 
the area of application calculations. These include a study 
of proton transfer in the HF trimer, which provides an analog 
of rearrangement reactions, and the extension of our previous 
work on Be and Mg clusters to Ca clusters. In addition, a 
very accurate investigation of the lowest-lying potential 
curves of the 0 2 molecule was completed. These curves are 
essential for evaluating different models of the terrestrial 
atmosphere nightglow. All these studies have been written up 
and accepted for publication (copies attached) . 

We have recently repeated some of our earlier 
studies of the Ne atom hyperpolarizability, stimulated partly 
by concern about very small uncertainties in the perturbed 
energies because of SCF convergence difficulties, and partly 
by the publication of a RASSCF investigation that predicted a 
very different hyperpolarizability. By calibrating our 
earlier calculations against full CCSDT results we can be 
confident that our earlier result (corrected for the 



convergence problems) should be much more reliable than the 
RASSCF result. This work has been written up and accepted 
for publication (copies attached) . 

Most effort this year has been devoted to a large- 
scale investigation of stationary points on the C 4 H 4 surface, 
and the thermochemistry and mechanisms of acetylene/acetylene 
reaction. Since acetylene is the last stable hydrocarbon 
species produced in the combustion of hydrocarbon fuels, it 
has been speculated that reaction between two C 2 H 2 species 
could be an important component of the overall combustion 
mechanism. Direct reaction of acetylene with itself to 
produce any C 4 H 4 products is likely to involve very high 
barriers; on the other hand, the very reactive vinylidene 
isomer of C 2 H 2 lies about 45 kcal/mol above acetylene. If 
significant quantities of vinylidene were formed these could 
react with acetylene to form a number of different C 4 H 4 
isomers. In fact, the existence of numerous minima on the 
C 4 H 4 surface, including some forms well-described only at the 
MCSCF level and thus not previously characterized, has 
seriously complicated the investigation. However, it seems 
likely that all relevant minima have been found, and some 
effort has been invested in locating saddle points on various 
pathways (if indeed there are barriers) . So far, it seems 
likely that vinylidene insertion into the acetylene CH bond 
will occur without a barrier, but there may be a barrier to 
CC bond insertion. This work is continuing. 
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Abstract 

Ab initio results are presented for the low-lying bound states of 0 2 that dis- 
sociate to ground-state atoms. Spectroscopic constants and dissociation energies 
are reported for the X 3 E a 1 A g , and b 1 E+ , c 1 E~, A/ 3 A„, M 3 E+ , 3 n u , 5 E U , 5 n u , 
and 5 n ff states. For the six lowest states, which have been experimentally charac- 
terized, we obtain accurate results at the multireference configuration interaction 
plus Davidson correction level. For example, we compute a D e value for the X 3 E ff 
state to be only 1.5 kcal/mole less than the experimental result after we have cor- 
rected for basis set superposition error. The 5 H ff state is estimated to have a D e 
of 0.16±0.03 eV, which suggests that the importance of this state in the nightglow 
should probably be reconsidered. 

I. Introduction 

The six strongly bound states of 0 2 that originate from ground-state atoms 
have been well characterized experimentally. More weakly bound states, such as 
5 n fl , have not been observed and will be difficult to characterize experimentally. Yet 
a knowledge of the potential energy curves of these states is important for modelling 
a number of phenomena. For example, the potential energy curves for these states 
are required for the calculation of transport properties and for modelling the energy 
flow in the Earth’s atmosphere. 

The terrestrial atmosphere nightglow is known to have two components: the 
f Mailing address: NASA Ames Research Center, Moffett Field, CA 94035 
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first is due to emission from the O 2 a 1 state, which is formed by the daytime 
photodissociation of O 3 , and the second is due to emission from a series of excited 
states, which are formed by the recombination of 0 atoms. This second process 
is very complex and is not completely understood. The first step involves the 
formation of excited O 2 molecules, 

0( 3 P) + 0 ( 3 P) + M -> O 2 + M. (1) 

These excited molecules can radiate, react, relax colli sionally, or dissociate. There 
seems to be general agreement that emission from several states, most notably 
is more intense than can be explained solely by a direct population mechanism and 
that a precursor must be involved. Both collision 

o; + M -> or + M (2) 

and reaction 

o; + o-+ or + o* (3) 

are known to be important energy transfer mechanisms for O 2 . For example, the 
reaction mechanism of eq. (3) is involved in the formation of 0( 1 5), which is re- 
sponsible for the nightglow green line. 

Wraight [x] and Smith [x] calculated the atomic association rate (i.e., the rate 
of eq. ( 1 )) into each of the bound states of O 2 and suggested that at 200 K 70% 
of all associations took place into the 5 n p state. The small D e of the 5 U g state 
leads to an increase in the dissociation rate state at higher temperatures, which 
results in a rapid fall off in O atom association rate at higher temperatures; their 
computed results were in good agreement with the experiments of Campbell and 
Gary. Hence Wraight suggested that the 5 II 5 state might play a very important role 
as a precursor in the formation of O 2 and 0( 1 iS') in the Earth’s atmosphere. 

To a large degree these conclusions were based on the 5 n fl potential energy curve 
calculated by Saxon and Liu, which has a D e of 0.22 eV. Unfortunately, the first- 
order configuration-interaction (FOCI) method used by Saxon and Liu yields only 
qualitatively reliable binding energies, while the calculated association rates are 
quite sensitive to the D e since the 5 H g lifetime must be sufficiently long that it can 
be quenched to a more tightly bound state. 
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Using a smaller estimate for the D e , Bates concluded that S II 9 was not an 
important precursor for the formation of the more strongly bound states or in the 
production of 0( 1 5). However, Bates suggested that it might play an important 
role in the formation of O3 in the lower atmosphere by the reaction 

0 2 * + 0 2 -+ 0 3 + o. (4) 

As these uncertainties in understanding the upper atmosphere are derived in 
part from the uncertainty in the dissociation energy of the 5 H 5 state, we have 
undertaken a study to compute this quantity very accurately. In addition, there is 
little information about the other weakly bound states of O2, such as 5 E U , which 
has been suggested as the state responsible for the perturbations in the A state. 
In this work we report on calculations on a number of the weakly bound valence 
states of O2. The calculation of binding energies, in particular, for such systems 
is a challenging problem, and we have attempted to estimate the reliability of our 
calculations by examining the convergence of our results for the six strongly bound 
states characterized accurately by experiment. The computational requirements 
for obtaining reliable spectroscopic constants for these states should provide useful 
information on the requirements for describing the hitherto uobserved states. We 
demonstrate that a multireference configuration-interaction (MRCI) treatment with 
a Davidson correction (denoted +Q) in an extended basis set yields reliable results. 

II. Method 

Most of the one-particle basis sets employed in this study are constructed 
using general contractions based on atomic natural orbitals (ANO) [1]. The first 
primitive set employed is that used in earlier work. The basis is derived from 
the (13s 8p) primitive set of van Duijneveldt supplemented with a (6d 4/ 2 g lh ) 
even-tempered set of polarization functions. The polarization functions are of the 
form o: = 2.5 n cto with n = 0...A: and 00 =0.13, 0.39, 1.24 and 2.61 for the 

d, /, g, and h functions, respectively. The geometric mean of the d exponents is 
derived from the optimum d exponent of Ahlrichs and Taylor. The primitive set 
was contracted to [6s 5 p 4 d 3/ 2g lh] based upon the natural orbitals obtained 
from atomic single and double excitation Cl calculations correlating both the 2s 
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and 2 p electrons, and employing an occupation number selection criterion of 10 -5 . 
Smaller sets, such as the [5s 4p 3d 2/ 1 g] set, are obtained by deleting the shell 
of orbitals with the smallest occupation number. A second contracted set was 
derived by supplementing the contracted set with a diffuse s (a = 0.076666) and p 
(a = 0.051556) function and uncontracting the outermost d function giving a set 
designated [5s 4p 3 + Id 2/ lflf]+(sp). Supplementing the basis with the diffuse 
functions improves the description of the oxygen atomic polarizability, an accurate 
description of which is necessary to calculate reliable interaction energies for weakly 
bound systems. 

To investigate the convergence with respect to the primitive basis set, we em- 
ployed an ANO set derived from a larger (18s 13p) primitive set. This set is quadru- 
ple zeta in the 2s space and yields an atomic SCF energy only 3 /xE# above the 
numerical Hartree-Fock result. This primitive set is supplemented with a (6d 5/ 4p) 
even-tempered set of polarization functions; cco = 0.13, 0.25, and 0.42 for the d, /, 
and g functions. The expanded / and g polarization sets are employed since work 
on Ne atom indicated that the smaller 4/ 2 g sets might not be saturated. The 
primitive set was contracted to [5s 5 p 3d 2/ \g\ based upon the natural orbitals 
from atomic calculations: natural orbitals with occupation numbers greater than 
10 -5 were selected. In addition, we employed the [9s 7 p 4 d 2/] segmented basis 
set from ref x. This set is a [7s 5 p] contraction of the (13s 8p) set of Partridge 
supplemented with two diffuse s, two diffuse p, and a (4d 2/) polarization set. The 
basis set was derived for a study of the high-spin states of NO and at the coupled- 
pair functional (CPF) level gives an O atom polarizability of 5.14 <Xq , or 95% of the 
recommended value. Only the pure spherical harmonic components of the d, /, g, 
and h functions are employed. 

The orbitals were optimized using the complete-active space self-consistent- 
field (CASSCF) approach. While the calculations are performed in D^h symmetry, 
f ull Dooh symmetry is imposed on the orbitals. In most cases the 2 p orbital and 
electrons are active. For the a 1 A g and b 1 E+ states, a state-averaged (SA) CASSCF 
approach is employed as they appear in the same D^h irreducible representation. 
MRCI calculations correlating both the 2s and 2 p electrons are performed: these 
include all single and double excitations from all configurations in CASSCF wave 
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function. The effect of higher excitations is estimated using a multireference analog 
of the Davidson correction, denoted +Q, and the size-extensive averaged coupled- 
pair functional (ACPF) method. Basis set superposition errors (BSSE) are esti- 
mated using the counterpoise method. All calculations were performed with the 
MOLECULE-SWEDEN program system. 

III. Results and Discussion 
A. Strongly bound states 

The calculated spectroscopic constants for the X 3 E“, a 1 A g and states 

in the different basis sets are compared with the experimental results in Table I. 
The r c and co e constants are derived from a fit in 1/r to the theoretical results 
on a 0.1 Uo grid. The D e values are computed using a supermolecule to mini- 
mize size-consistency errors. Overall the results are in rather good agreement with 
experiment. For the ground state expanding the basis set from [4 s 3 p 2d 1/] to 
[55 4 p 3 d 2/ lg] results in little change in r e and u> e , but D e is increased by almost 
0.2 eV. Expansion to [6s 5 p 4 d 3/ 2 g 1 h] has a much smaller effect, increasing D e 
by about 0.05 eV. These results suggest that the calculated spectroscopic constants 
in our biggest basis set are near the basis set limit. Extending the [5s 4 p 3 d 2/ lg] 
basis set with diffuse functions (and uncontracting the outermost d primitive) in- 
creases the calculated D e by about 0.025 eV, which is considerably more than the 
0.009 eV obtained by substantially expanding the primitive basis set. Based on 
these results we use the [5s 4p 3 + Id 2/ lg] + (sp) basis as our standard basis 
set in the remainder of this work. In this set the r e and u? c values are in excellent 
agreement with experiment. At the MRCI level the calculated D e values are about 
0.15 eV less than the experimental results, whereas including the Davidson correc- 
tion for higher excitations gives D e values about 0.03 eV less than the experiment. 
This differs from previous work employing a [5s 4 p 3 d 2/ lg] ANO basis set, which 
suggested that the MRCI+Q level of treatment in these basis sets overestimated 
D e . The discrepancy occurs since the highest a g and a u orbitals were deleted in the 
earlier work. At long range this corresponds to deleting a combination of the fourth 
and fifth s ANOs, while near r e the orbitals are a combination of all possible a 
components. Deleting these orbitals artificially raised the atomic asymptote (when 
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twelve electrons are correlated) resulting in a 0.08 eV overestimation of D c . In com- 
parisons with double zeta plus polarization basis eight-electron full Cl calculations, 
the MRCI+Q treatment was shown to overestimate D e by about as much as MRCI 
underestimates it. However, when 2s correlation is included it is unlikely that the 
MRCI+Q treatment will significantly overestimate D e and we consider the results 
at this level as the most reliable. In the ANO basis that includes h functions the 
MRCI+Q D e value is slightly larger than experiment, but the difference is of the 
same order as the superposition error, as discussed below. The reliability of the 
Davidson correction is corroborated by the very similar results that are obtained 
at the ACPF level. Indeed, the ACPF results would suggest that the Davidson 
correction slightly underestimates the effect of higher excitations. 

Compared with the X 3 £“ MRCI results with eight electrons correlated of Ref. 
x, the effect of 2s correlation is to reduce D e by about 0.1 eV. This is considerably 
larger than the 0.04 eV reported earlier; as noted above the error in the older 
twelve-electrons correlated results arose from truncating the virtual orbital space. 
However, even this larger value is substantially less than the 0.3 eV reduction for 
N 2 . The reduction in D e in that case has been attributed to single excitations in 
the N atom from 2s to 3d, with a concomitant spin and symmetry recoupling of 
the 2 p electrons. In N 2 this effect disappears at r e , because all the 2 p electrons are 
coupled to a singlet to form the triple bond. The same atomic correlation effect 
arises in the O atom, but in the O 2 ground state part of the effect persists even 
at r e , because the 2p electrons are not completely recoupled in forming the bond. 
Hence the reduction in D e from 2s correlation will be less in O 2 than in N 2 - 

The basis set superposition errors, BSSE, for the different basis sets were esti- 
mated using the counterpoise method. At the MRCI+Q level the BSSE (per atom) 
at r=2.3 a 0 is 0.022 eV for the [5s 4p 3d 2/ 1 g] and [5s 4p 3 + Id 2/ lg]+(sp) 
basis set and 0.019 eV for the [5s 5 p 3d 2/ 1^] basis set. Most of the reduction in 
the BSSE is due to the improvement of the primitive polarization basis; expanding 
the s and p primitive spaces reduced the BSSE very slightly. Thus, corrected for 
BSSE the calculated D e of the X 3 E~ state is 0.063 eV (1.45 kcal/mole) smaller 
than the experimental D e . We estimate that about half this error is due to basis 
set deficiencies, based on our experience with N 2 . Core (Is) correlation is likely to 
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contribute on the order of 0.7 kcal/mole. 

In Table II, we compare the spectroscopic constants at the MRCI+Q level in 
the [ 5 s 4 p 3 + Id 2 / l#]+(«sp) basis set for the six bound states, dissociating to 
ground state atoms, that have been experimentally characterized. The calculated 
results are in excellent agreement for all of the states. Except for the a 1 A g and 
b states the calculated D e are all within 0.025 eV of the experimental results. 
For the a 1 A g and b states the error is 0.053 and 0.041 eV, respectively. We 
have investigated the possibility that the use of SA-CASSCF orbitals is responsible 
for this, but the error from this source is only 0.01 eV. Overall, our results indicate 
that the computational procedure and basis set employed gives a good description 
of all of the low-lying states. In fact, the calculated potentials for the c 1 E“, A ' 3 A u , 
and A 3 E^ states, tabulated in Table III and plotted in Fig. 1, agree well with the 
RKR potentials, indicating that we are obtaining a good description for all r values. 


B. Weakly bound states 

The computational requirements necessary to obtain an accurate description 
for weakly bound systems are severe, because the errors associated with the one- and 
n-particle space truncation can be of the same order of magnitude as the well depth. 
It is important that the computational procedure be able to describe the atomic 
polarizabilities well and to minimize the BSSE. The ANO basis sets employed in 
this work were contracted based upon a correlated treatment of the 0( 3 P) state. 
Thus, it is expected that, in some sense, these basis sets minimize the BSSE. In 
addition, by uncontracting or supplementing with diffuse functions, the basis sets 
permit an accurate description of the atomic properties. 

For the weakly bound states of O 2 dissociating to ground state 0( 3 P) atoms, 
the long range form of the potential is of the form 


-V LR (r) - — + — + — + 


where the coefficient C 5 arises from the quadrupole-quadrupole interactions, C 6 is 
the dispersion coefficient arising from the induced dipole-induced dipole interaction, 
and C& describes both dispersion and induction. We note parenthetically that the 
free O atom does not have a quadupole moment, of course: it is only when the 
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atoms are perturbed, such as by another atom, that the degeneracy is lifted and it 
is meaningful to speak of the quadrupole moments of a specific configuration. In 
recent work on the high-spin states of NO we examined whether it is important 
to constrain the 0 atom at long range to be spherical, i.e., to ensure that the all 
of the states dissociating to ground state atoms have the same asymptotic energy. 
We found that these constraints were not necessary for treatments that included 
correlation. In addition to substantially increasing the length of the Cl expansion, 
the constraints reduced the relibility of the Davidson correction in the vicinity of 
the minimum. Thus the constraints resulted in larger calculations that did not 
improve the description of the long-range interaction, and degraded the description 
of the short-range interaction. We do not require the 0 atoms to be spherical at 
long range here, but use a supermolecule to define the asymptotic energy for each 
state. 

C 5 is positive for the ^g, 3 II U and 5 IIg states. Hence at long range this is the 
dominant term and these states are noticably more bound than the other states of 
O 2 — see Figs. 1 and 2. In fact, for r values longer than about 5 ao, the 5 n s state 
is the lowest-lying state of 02- This, coupled with the high degeneracy of the state, 
led to the speculation concerning the importance of this state as a precursor for the 
formation of the state and 0( 1 5) states. As discussed by Wright and in detail 

by Bates, the importance of this state as an intermediate is strongly dependent upon 
its D e value. Using the [5s 4p 3 + Id 2/ lp]+(sp) basis set employed for the other 
states of 0 2 , we obtain a D e of 0.131 eV at the MRCI+Q level. This is significantly 
less than the 0.22 eV obtained at the first-order Cl (FOCI) level. Because of the 
possible significance of this state in explaining the details of a number of processes, 
it is important to calibrate the accuracy of our calculations. 

The results of our calibration study for the 5 n s state are given in Table IV. 
First, we note that the inclusion of 2s correlation substantially deepens the well 
and decreases r e . The difference with the effect on the ground state occurs be- 
cause the atomic correlation terms involving 2s to 3 d excitation with a recoupling 
of the 2 p electrons are still present at r e . In fact the 2p recoupling allows some 
additional bonding, thus 2s correlation increases the bonding energy for the 5 1I 9 
state. Second, while the magnitude of the Davidson correction is similar to that for 
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the experimentally characterized states, the relative magnitude of the correction is 
substantial — it more than doubles the well depth. While corrections of this type 
and of a similar magnitude have been shown to be accurate for other weakly bound 
systems, it is important to calibrate the correction for this state, as full Cl stud- 
ies have demonstrated that the Davidson correction can overestimate the effect of 
higher excitations. This was done by performing expanded reference calculations at 
the eight-electron level using the [45 3 p 2d 1/] basis set. A second p shell was added 
to the CASSCF active space, CASSCF(2p,2p f ), and then an MRCI was performed 
with a reference list consisting of all CSFs that could be formed by distributing 
the eight electrons in all ways among the 3<x gy 4 cr g , 3 <r u , l7r u , 2ir u and liz g orbitals. 
That is, the antibonding orbitals of the second p shell are not included in con- 
structing the reference space. The MRCI based on this expanded reference space 
is considerably (0.021 eV) deeper than the smaller MRCI, but the MRCI+Q result 
is only 0.002 eV shallower than the previous MRCI+Q result. In addition, using 
the smaller active space and the [5 s 4p 3 + Id 2/ lg]-\-(sp) basis set, the ACPF and 
MRCI+Q results are again nearly the same — the ACPF D e is 0.003 eV smaller. 
These results suggest that even at the eight-electron level the MRCI+Q results are 
reliable for the b IL g state. Since the importance of higher excitations will increase 
when 25 correlation is included, it is likely that the +Q correction is actually an 
underestimate at the twelve-electron level. This is supported by the ACPF results, 
which yield a D e 0.021 eV larger than the MRCI+Q result. If the 2s orbitals and 
electrons are included in the active space, CASSCF(25,2p), and this CASSCF used 
as the reference space, then /^increases by 0.011 and 0.015 eV at the MRCI+Q 
and ACPF levels, respectively. However, based on accurate calculations on N 2 , we 
believe that the CASSCF(2.s,2p) reference MRCI+Q and ACPF calculations may 
overestimate the effect of higher excitations somewhat. Nevertheless, we can use 
these results to obtain an error estimate for the n-particle treatment by assuming 
that the CASSCF(2p)/MRCI+Q is a lower bound and the CASSCF(25,2p)/ACPF 
is an upper bound. Correcting for BSSE we obtain a D e of 0.139 ± 0.018. Basis set 
improvements will unquestionably increase the calculated D e . For the X 3 E 5 state, 
by comparing our best computed result to experiment we observe that the basis 
set incompleteness is about twice the BSSE, so it could be argues that rather than 
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subtracting the BSSE we should add it to our computed D e , to account for basis 
set incompleteness. Since this procedure is far from rigorous we have increased our 
uncertainty estimate by 0.01 eV to account for the uncertainty in our one-particle 
extrapolation. This leads to our best estimate of 0.16iL0.03 eV for the D e of the 5 Hj 
state. This value is slightly larger than our preliminary estimate, quoted by Bates, 
of 0.14±0.03 eV. While our estimate is smaller than the FOCI result employed by 
Wraight, it is large enough that it might be necessary to reconsider the importance 
of the 5 U g state as a precursor. 

In Table IV we compare our computed spectroscopic constants with the FOCI 
results of Saxon and Liu. Our calculated r e is 0.2 ao longer and u> e is much smaller 
than the FOCI results — see Fig. 3. The difference is a consequence of the FOCI 
procedure, not the one-particle basis set or method of optimizing the molecular or- 
bitals, as we obtain similar FOCI results using our basis sets and CASSCF orbitals. 
The FOCI approach substantially overestimates the interaction energy in the region 
of the potential well: the curve is not even qualitatively correct. In fact, the FOCI 
approach overestimates the interaction of many of the weakly bound states in this 
region. For NO + we also found states for which the FOCI binding energies were 
substantial overestimates. 

The spectroscopic constants for the other weakly bound systems are summa- 
rized in Table V. The calculated spectroscopic constants for the 5 E“ and 5 II U states 
are in surprisingly good agreement with the FOCI results. However, this agreement 
must be fortuitous, because all of the binding in these van der Waals bound states 
arises from the dispersion interaction and the FOCI does not include the simul- 
tanous single excitations required to describe it. Based on previous work on the 
high spin states of N 2 and NO, we estimate that our calculated D e values are accu- 
rate to within 20%, with the true potential being deeper than the calculated result. 
Because of the positive C 5 term, the 3 H U and ^ states have slightly larger D e 
and shorter r e values than the weakly bound states, where C$ is the first non-zero 
term in the multipole expansion. 

Borrell et al. found perturbations in the N= 9, 11, and 13 rotational lev- 
els of u=ll for the A 3 £+ state. These rotational levels lie between 61±15 and 
15±15 cm -1 below the dissociation limit. They suggested, based on the poten- 
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tials of Saxon and Liu, that the 5 S“ state was the most likely candidate for the 
perturbing state. Our ACPF D e value of 47 cm -1 , which is expected to be up to 
20% too small, is not inconsistent with this assignment. We note that the 1 II S and 
3 n„ states cross the A 3 £+ state in this region as well. However, if they were the 
perturbing states, they have would have been expected to affect levels below N— 9, 
since they have significantly larger D e values than the 5 n s state. Thus our results 
agree with the hypothesis of Borrell et al. that 5 £ ” is the most probable candidate 
for the perturbing state. 

IV. Conclusions 

We have demonstrated that a CASSCF(2p)/MRCI+Q description in a 
[5^ 4p 3d 2/ lg] ANO basis set supplemented with diffuse functions provides a 
quantitative description of the six lowest states of 02- The calculated potentials 
are within 0.05 eV (1.5 kcal/mol) of the (accurate) experimental results. In addi- 
tion, we have investigated the importance of substantially expanding the primitive 
basis set and demonstrated that such expansions yield insignificant improvement in 
the spectroscopic constants: the expanded basis set recovers 0.57 eV more atomic 
correlation but D e increases by only only 0.01 eV. 

Potential energy curves have also been reported for the weakly bound states 
of O 2 . Based on a series of calibration calculations the 5 H S state is estimated to 
have a D e of 0.16±0.03 eV. This estimate is slightly larger than our preliminary 
estimate, quoted by Bates, but the upper bound of our D e is sufficiently large that 
the importance of this state as a precursor should be reconsidered. The computed 
spectroscopic constants for the other weakly bound systems are in reasonably good 
agreement with previous FOCI estimates. 
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Table III . Calculated energy values for O2 (in cm *). 


r 


A ' 3 A u 


5 n 

Ai s 

3 n u 

% 

5 s- 

5 n u 

2.10 

25784.41 

27954.91 

29276.27 






2.30 

5491.75 

7408.75 

8493.08 






2.50 

- 4248.81 

- 2452.53 

- 1560.47 


21041.97 




2.70 

- 8143.70 

- 6374.92 

- 5638.27 






2.80 

- 8793.02 

- 7011.36 

- 6340.75 

15902.33 

8111.17 




2.90 

- 8893.33 

- 7085.20 

- 6474.02 






3.00 

- 8615.03 

- 6770.91 

- 6212.59 



30511.30 



3.01 


- 6723.72 

- 6170.39 

7334.67 

4240.44 




3.10 

- 8089.77 

- 6206.20 

- 5694.97 






3.30 

- 6676.75 

- 4729.16 

- 4296.81 






3.50 

- 5187.38 

- 3250.95 

- 2883.00 

- 45.36 

2016.92 

9530.66 


6620.31 

3.70 

- 3881.34 

- 2066.34 

- 1754.19 

- 778.25 





3.80 




- 943.19 





3.90 




- 1025.92 





4.00 

- 2424.63 

- 985.14 

- 749.31 

- 1052.98 

1169.42 

3444.33 

4171.39 

2342.29 

4.10 




- 1043.47 





4.30 




- 964.44 





4.40 







1752.90 


4.50 

- 1092.02 

- 343.33 

- 213.53 

- 853.44 

297.38 

1048.55 


752.11 

4.70 




- 739.23 





4.80 







694.87 


5.00 

- 510.45 

- 170.19 

- 94.81 

- 584.91 

- 90.20 

186.92 


195.53 

5.50 

- 252.66 

- 99.45 

- 52.70 

- 388.95 

- 188.69 

- 80.91 

92.17 

19.57 

5.75 





- 190.53 

- 122.98 


- 11.48 

6.00 

- 133.66 

- 61.05 

- 30.06 

- 258.31 

- 179.19 

- 136.78 

- 5.68 

- 25.15 

6.25 






- 135.23 

- 20.94 

- 29.62 

6.50 

- 75.86 

- 38.09 

- 16.58 

- 173.11 

- 142.72 

- 125.97 

- 26.51 

- 29.23 

6.75 







- 27.07 


7.00 






- 100.09 

- 24.95 

- 23.48 

7.50 

- 29.88 

- 16.16 

- 5.45 

- 81.80 

- 78.28 

- 75.64 

- 19.61 

- 16.96 

8.00 







- 13.99 

- 11.79 

8.50 

- 15.02 

- 8.91 

- 3.00 


- 43.00 


- 9.67 


9.00 


- 7.17 

- 2.78 





- 7.25 

10.00 

- 7.30 

- 5.14 

- 2.66 

- 19.03 

- 19.94 

- 19.83 

- 3.49 

- 3.06 

12.00 







- 2.88 


20.00 







0.00 

0.00 

50.00 


0.00 

0.00 

0.00 

0.00 

0.00 




12 


Table IV. Calibration study for the O 2 5 II fl Spectroscopic Constants 0 . 


05 

CM 


CM 


^ O « O 

' 05 s — ' 

<n t— 1— I 05 05 CO 

CO to o O VO 

<u O r— i * i—i O 

C) o o 3^ o © 


CO 


o ^ 


CM 


t- 

CN 


o co 00 _ . ^rr' 

COOOcOlC 1 ^^^^ 

H ^ HH lOCOCOlO CM 
• ©••rHOrH© CM 
o • • • • • * 


£ 
6 CM 

3 CM 


o o 

CO CO 
s CM CM 

TfH ✓ 05 s " 

05 t- 05 CM CM 

hO^COcm^ 

' — * t—h t— h ' — CM 


CM 

CM 

O 

t— 

CM 

CM 


t- 

co 

CM 


t-H 


CO 

CO 


CO 

0 


0 


^ — N 

Tf ^ 

N ' ' 

to _ 

^ CO 

0 O 

CO 0 

^ CM 

d co 

rH O 

rH 0 

CM 

* T— 1 

co ■ 

y • 

• 

• 


^ Tf 

Tf ^ 


r~- 

05 

CO 


fa 
fa 
O 

'fa 'fa 
CM^ CM^ 

«r «r 

CM CM 

fa fa 

. r o o 

ce^ co co 

_? fa CO CO 

^ O <1 <! 

H H <| O O 
fa fa | | | 

CN AAA 
fa CO «0 fa fa fa 
CO CM CM «5 «5 «5 


fa 

fa 

O 

< 


fa 

fa 

o 

c 

I 


+ + + + + 
O 5 O 5 fa fa 


^7 e*. * + ~» * 4 ~ , » * + ~* H -> 

®><N M N M M 


+ 

1 1 

H-*' 

CM 


O 

tO 

CO 


o 

oo 

co 


o 

u 

0 

15 

13 

1 

CM 


^ ^ ^ ^ ^ ^ 

CM CM 1 — It — l CM CO 

"d ’"d *d ’"d "d 

^ CM CM CO 


"d *d *d *d 'd 


+ + + + + 

CO co CO CO CO 


+ 


fafafafafafafafafafafa 

NNCOCO^lO^^^^T^ 


co co 
05 05 


coco«0«5«5CO«0«0 

TfiOOlOtOtOlOlO 


CM 

^d 

CO CO 

fa fa 

^ to 

«0 co 

to to 


0 

fa 

•d 

d 

ce 

d 

o 

£ 

CO 


00 co co d 

CO CO 05 00 

o o o o 

© o © © 

00 05 H CO CO 

Tf t- CO 05 to 

O © © © © 

00000 


^ ^ 10 

to CO 00 


to 

00 


to 00 CO to co 
o co CO 00 co 


CM 00 to d 

CO 05 d 

CM CM © rH 

d ^ 

H C5 00 

00 to d o co 

TjJ CM CO rH CO 

Ttf d Tf Tj5 


o 

'—pH 

fa s 

cm y 

fa 

C*" y 

fa I 

O + 

CO 


fa 

fa 

O 

fa 

«5 

+ 


CO 

< 

O 

^ JU 

*+-> ^ 


<55 fa 


CM 

*d 


§ ^ + 
M m « 

u ft. 

a; co co Tf 

15 «o <0 <0 

00 ■^.. T ^.. LO 


^ > 4 ^ 

^ "d 

CO CO 

fa fa 
^ to 

co <0 

to to 


fa 

CM 

fa 

O 

CO 

CO 

< 

O 

D 

fa 

-+■* 

bO 

.5 

* C/2 

d 

w 

15 

w 

15 

fa 

h-^ 

a 

15 

5h 

<3 

fa 


d 

0) 

> 

’So 

w 

r-H 

d 

w 

u 

a 

+ 

hH 

o 


05 

fa 


w 

3 

CO 

0) 

(H 

t— I 

o 

PH 


15 

fa 

H-3 

05 

M 

03 

D 

w 

05 

fa 

"d 

05 

o 

a 

05 

u 

15 

fa 


> 

15 

05 

I>- 

O 

o 


O' 

+ 

0 S 

5 a 

11 a? 

s- fa 
bO ^ 
d r£] 

*2 ^ 
p CO 

"d 

15 


ffl 


15 
O 
0 
15 

,15 

fa 13 

8 * 
W «3 
W w 
a c3 


ci 

3 

o 

13 

15 


"d 

15 

-4-Z 

u 

15 

M 

o 

O 



d The MRCI binding energy at r=4.0 a 0 is 0.054 eV. 
e Corrected for BSSE, the MRCI+Q D e is 0.121 eV. 
f Corrected for BSSE, the MRCI+Q D e is 0.122 eV. 


lO 



Table V. Spectroscopic constants for the weakly bound states of O 2 . 


Calculation a 

r e (a 0 ) 

We (cm J ) 

D e (cm *) 

5 n u 

[5s 4p 3 + Id 2/ lg] + ( sp ) 

[5s 4p 3 + Id 2/ 1$] + (sp)-ACPF 

6.71(6.63) 

6.66 


36.8(37.3) 

36.7 

Saxon-Liu 6 

7.2 


34.0 

5 s- 

[5s 4p 3 + Id 2/ 1 g\ + (sp) 

[5s 4p 3 + Id 2/ lg] + (sp)-ACPF 

6.49(6.24) 

6.15 


31.1(42.8) 

46.8 

Saxon-Liu 

6.5 


39.5 

3 n u 

[5s 4p 3 + Id 2/ 1 5 ] + (sp) 

5.83(5.65) 

60(59) 

147.0(191.6) 

Saxon-Liu 

5.5 


165.2 

% 

[5s 4p 3 + Id 2/ 1 g] + (sp) 

6.17(6.09) 

57(58) 

111.5(137.8) 

Saxon-Liu 

6.2 


117.9 


“The numbers in parenthesis include a Davidson correction. 
6 Ref. xx 
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Introduction 


There has been considerable recent interest in the properties of small clusters (see 
for example Refs. 1-10), motivated principally by two issues. The first is the question 
of convergence of cluster properties towards the bulk values. Of course some properties 
will approach the bulk value more quickly than others as the cluster size is increased. The 
second issue is interaction between theory and experiment. The study of small clusters has 
progressed very rapidly since accurate experimental studies may be used to evaluate the 
predictive reliability of different theoretical methods, and then accurate theoretical studies 
may be used to evaluate or aid in the design of new experimental techniques. 

Our studies of small clusters 11 ’ 12 have focused on computing the structures, binding 
energies, vibrational frequencies and infrared intensities of the trimers and tetramers of the 
alkaline-earth elements beryllium and magnesium using elaborate treatments of electron 
correlation. We have also examined 13 the equilibrium structures and binding energies of 
the pentamers, Be 5 and Mg 5 . In the present study we extend our investigations to in- 
clude Ca 3 and Ca 4 . There are three previous studies 14 ’ 15 of Ca 4 (and none of Ca 3 ) that 
have incorporated electron correlation effects. In one Ca 4 study a total binding energy of 
18.3 kcal/mol, relative to four Ca atoms, was obtained at the single-reference single and 
double excitation configuration interaction (CISD) level of theory, including Davidson’s 
correction 16 for higher excitations. In the other, a binding energy of xx.x kcal/mol was ob- 
tained using a multireference Cl approach (based on SCF orbitals). Based on our previous 
studies of the Be and Mg clusters, it is likely that even the higher value is a substantial 
underestimate of the true binding energy. 

Another important component of our earlier studies 11-13 ’ 17 is the comparison of ge- 
ometries and binding energies obtained with various electron correlation methods. The 
s-p near-degeneracy effects in the alkaline-earth valence shell are very large, and strongly 
influence the binding in the small clusters. It is not only essential to describe these non- 
dynamical effects accurately, however, but also to account properly for dynamical correla- 
tion in order to obtain reliable binding energies for these systems. Hence the most desirable 
treatment might appear to be a full valence complete-active-space SCF (CASSCF) calcu 
lation followed by a multireference configuration-interaction (MRCI) calculation. This is 
indeed an excellent level of treatment, but unfortunately it becomes very expensive to 
apply in large basis sets and at many geometries. 
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In our studies of the lighter alkaline-earth clusters, we have made extensive use 
of the single and double excitation coupled-cluster approach (CCSD) corrected with a 
perturbational estimate of triple excitations (CCSD(T)). 18 The CCSD(T) method performs 
very well in comparison with MRCI results for the lighter alkaline-earth clusters, and 
appears to treat both the non-dynamical and dynamical correlation effects in these systems 
accurately. We have also established the reliability of the CCSD(T) method by full Cl 
comparisons on Be 3 . 20 compared results obtained at the MRCI level of treatment with 
those obtained using the single and double excitation coupled-cluster (CCSD) approach 
and the CCSD(T) method, in which a perturbational estimate of the effects of connected 
triple excitations has been included. 18 In these studies we established that the CCSD(T) 
method performs very well in reproducing results obtained from the MRCI approach. 
Based on many comparisons of full configuration interaction (FCI) and MRCI energies, 19 
it has been asserted that the MRCI results should be very close to the n-particle limit, 
and we have established 20 that the MRCI and CCSD(T) results for Be 3 are indeed very 
close to the analogous FCI quantities. For comparison purposes, we have again used the 
CCSD, CCSD(T) and MRCI methods in examining the potential energy surfaces of Ca 3 
and Ca 4 so that it may be determined how well the coupled-cluster methods perform for 
the series of alkaline-earth metals. 

In the next section we describe the computational methods employed in this study 
and in the following section our results are presented and discussed. A comparison of our 
new results for the Ca clusters and our previous results for the Be and Mg clusters is also 
given. The final section contains our conclusions. 

Computational methods 

Two atomic natural orbital 21 (ANO) basis sets have been used in this study. The 
(22s 17p) primitive basis set is that of Partridge 22 and was augmented with a (4 d 3/) even 
tempered polarization set defined by a = 2.5 n a„ for n = 0,...,k. The a 0 values for the d 
and / functions are 0.0232 and 0.0440, respectively. The smallest basis set consists of 5s, 

4p and Id ANOs and will be denoted [5s Ap Id}. The larger basis consists of 6s, 5p, 2d 
and If ANOs and will be denoted [6s 5p 2d If}. The ANO contraction coefficients were 
obtained by averaging the natural orbitals from single and double configuration interaction 
(CISD) calculations on the lowest and 1 P states of atomic Ca. Only the pure spherical 
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harmonic components of the d and / functions have been used. 

As discussed in the Introduction, the CCSD, CCSD(T) and MRCI methods have 
been used to treat electron correlation. In all cases, only the Ca 4s electrons have been 
included in the correlation procedure. The coupled-cluster wave functions are based on 
self-consistent field (SCF) molecular orbitals while the MRCI wave functions are based 
on CASSCF molecular orbitals. All valence electrons were allowed variable occupancy in 
all valence orbitals in the CASSCF calculations. References for the MRCI wave functions 
were selected using a 0.05 threshold — that is, all occupations having a component spin- 
coupling with a coefficient of 0.05 or larger in the CASSCF wave function were used as 
references in the MRCI procedure. 

In analogy with small Be and Mg clusters, the equilibrium geometries of Ca 3 and 
Ca 4 were constrained to have D 3 h and Tj, symmetry, respectively. Harmonic frequency 
analyses demonstrate explicitly that these geometries are indeed minima on the Ca 3 and 
Ca 4 potential energy surfaces. In addition, a linear structure for Ca 3 was optimized and 
found to be significantly higher in energy than the equilateral triangle. Hence it is expected 
that the equilateral triangle and tetrahedron are the global minima on the Ca 3 and Ca 4 
potential energy surfaces, respectively. 

The quadratic, cubic and quartic force constants of Ca 3 and Ca 4 have been deter- 
mined numerically and are given in symmetry internal coordinates. The symmetry internal 
coordinate definitions are: 

Ca 3 


Ca 4 


SM) = 

5 a .(e') = 
5 2 *(e') = 

Si(ai) = 

S 2a (e) = 
S 2b (e) = 


—/=( r l + r 2 + 7 * 3 ) 
v3 


Ve 


(2 r a - r 2 - r 3 ) 




v5 (r 1 +r J +r 3 +r 4+ r 5 +r,) 


Vl2 

1 


(2 7*1 - r 2 - r 3 + 2r 4 - r 5 - r 6 ) 


V4 


(r 2 - r 3 + r 5 - r 6 ) 


( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 
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( 7 ) 


S 3l (< 2) = ^=(r 2 -r 5 ) 

S 3y (t 2 ) = -^(r 3 -r 6 ) (8) 

Ss*(*a) = ^=(ri-r 4 ). 

For Ca4, the numbering of bonds is such that if r 3 connects one pair of atoms, 7-4 connects 
the other pair, and similarly for the pairs {72,7-5} and {7-3, re}. According to convention, 
the a component of the doubly degenerate coordinates is defined such that it is symmetric 
with respect to a cr v reflection plane, whereas the b component is antisymmetric. The 
x, y and z components of the triply degenerate coordinate are also defined according to 
established convention. 24 

The precision of the central difference numerical procedures used to obtain the force 
constants has been closely monitored: the uncertainty in the harmonic frequencies should 
be less than 0.1 cm -1 and in the fundamental frequencies less than 0.5 cm -1 . The anhar- 
monic analyses have been performed with the SPECTRO program package which uses 
second-order perturbation theory; Ca 3 has been treated as a symmetric top 23 and Ca4 has 
been treated as a spherical top. 26 The coupled-cluster calculations were performed with the 
TITAN set of programs 27 interfaced to the MOLECULE-SWEDEN suite of programs. 
The MRCI calculations were performed with MOLECULE-SWEDEN. 

Results and Discussion 

A. Equilibrium structures and binding energies 

Table 1 contains the equilibrium bond distances, rotational constants and dissociation 
energies (atomization energies) of Ca 3 and Ca4 computed in this study. For Ca 3 the 
CCSD level of theory substantially underestimates the D e value (i.e., 5.7 kcal/mol or 47% 
compared to MRCI) and therefore yields a bond length 0.22 ao too long. However, the 
CCSD(T) level of theory exhibits a significant improvement over the CCSD results: D e 
is only 0.9 kcal/mol less than the MRCI result and the equilibrium bond distance differs 
from the MRCI value by only 0.01 a 0 . A similar situation is found with Ca 4 — the CCSD 
level of theory substantially underestimates the binding energy and gives an equilibrium 
bond distance 0.11 a 0 too long. However, as was found for the Be and Mg clusters, 11 - 12 
the CCSD level of theory performs better for Ca 4 than it does for Ca 3 . The CCSD(T) 
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results for Ca 4 (using the smaller ANO basis set) are in very good agreement with the 
respective MRCI quantities. The CCSD(T) D e is only 0.2 kcal/mol less than the MRCI 
value and the CCSD(T) r e is actually 0.02 a 0 shorter than the MRCI r e . As with our 
earlier studies 11 - 12 of small Be and Mg clusters, it thus appears that the CCSD(T) level of 
theory provides a very good description of electron correlation effects in small Ca clusters. 

It is interesting to note that comparison of the equilibrium bond distances obtained 
with the CCSD, CCSD(T) and MRCI methods suggests that the bonding in small Ca 
clusters is intermediate between the bonding in small Be and small Mg clusters, as would 
be expected based on the bulk binding energies. 29 For Be 3 , where sp hybridization is known 
to play an important role in the binding, the CCSD level of theory yields a reasonable 
equilibrium bond distance when compared to MRCI or CCSD(T). However, for Mg 3 , where 
the binding is more dominated by dispersion, the CCSD level of theory gives a bond 
distance that is significantly too long (0.55 ao). The discrepancy in the CCSD bond length 
for Ca 3 is between that found for Be 3 and Mg 3 , but is much closer to the value obtained 
for Be 3 than that found for Mg 3 . It therefore seems that sp hybridization is an important 
component of the bonding in Ca 3 , though not as important as it is in the bonding of Be 3 . 

For Ca 4 , only the CCSD and CCSD(T) levels of theory could be used in conjunction 
with the larger ANO basis set since the MRCI procedure would have been prohibitively 
expensive. Using the [6s 5 p 2d If] ANO basis set, the CCSD(T) equilibrium bond distance 
for Ca 4 is 0.29 ao shorter than the analogous Ca 3 value indicating the increased importance 
of sp hybridization in the bonding as the cluster size becomes larger. The best D e values for 
Ca 3 and Ca 4 obtained in this work are 12.1 kcal/mol and 31.5 kcal/mol, respectively. As 
expected, our best computed D e for Ca 4 is substantially larger than the previously best 
computed value. 14 Based on the fact that the separated-atom limit is described better 
than the molecule, the D e values for Ca 3 and Ca 4 are probably underestimated somewhat. 
We say “probably” because it is not certain whether the effects of core-correlation will 
increase or decrease the D e values. However, based on a recent study of Ca 2 by Dyall and 
McLean, 30 it is likely that the effects of core-correlation will not affect the D e values by 
more than 1-2 kcal/mol. 


B. Vibrational frequencies of Ca3 

The quadratic force constants and harmonic frequencies of Ca 3 obtained at the 
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CCSD, CCSD(T) and MRCI levels of theory are presented in Table 2. The [6s 5 p 2d 1/] 
ANO basis set was used with all of these methods. In view of the underestimation of the 
bond strength at the CCSD level, it is not surprising that the CCSD quadratic force con- 
stants and harmonic frequencies are noticeably smaller than the analogous MRCI values. 
Conversely, the CCSD(T) quadratic force constants and harmonic frequencies are in excel- 
lent agreement with the MRCI quantities, being only 1 cm 1 and 2 cm smaller for the 
a[ and e f modes, respectively. The best computed harmonic frequencies of Ca 3 obtained 
in this work are 95 cm" 1 and 85 cm" 1 for the a\ and e' modes, respectively. These values 
should be reasonably reliable and are probably low rather than high, assuming that the 
binding energy is still somewhat underestimated. The harmonic frequencies of Ca 3 are 
small not only because the bonding in Ca 3 is relatively weak, but also because of the large 
mass of the Ca atom and the l/y/rn dependence of the harmonic frequency, where m is 
the reduced mass of the system. 

We have computed cubic and quartic force constants using the CCSD(T) method, 
since the CCSD(T) and MRCI harmonic frequencies are in excellent agreement and the 
former approach is significantly cheaper. The complete set of cubic and quartic force con- 
stants and the resulting anharmonic constants are presented in Table 3. Table 4 contains 
the fundamental vibrational frequencies of Ca 3 determined via second-order perturbation 
theory. In addition, Table 4 presents the infrared (IR) intensity of the e' vibration deter- 
mined using the double harmonic approximation. 

The absolute anharmonic contribution to the vibrational frequencies is small, only 
2 cm -1 for both vibrations, and the percentage effect is about half that observed previously 
for the Mg clusters. Combining the MRCI harmonic frequencies with the anharmonic 
correction obtained with the CCSD(T) method gives 93 cm 1 and 83 cm as our best 
estimates for the fundamental vibrational frequencies of the a\ and e' modes, respectively. 
Because the IR intensity of the e' mode is so small, the best prospect for determining the 
vibrational frequencies experimentally is likely to be an indirect technique such as negative 
ion photoelectron spectroscopy. 


C. Vibrational frequencies of Ca4 

The quadratic force constants and harmonic frequencies of Ca 4 , determined at the 
CCSD, CCSD(T) and MRCI levels of theory, are presented in Table 5. A noteworthy point 
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is that the CCSD harmonic frequencies are in better agreement with the MRCI values than 
was found for Ca 3 . This is especially true for u? 2 and u; 3 where the differences are only 
4 cm -1 and 8 cm -1 , respectively. This observation supports our earlier analysis regarding 
the increased importance of sp hybridization, and consequently covalent bonding, in Ca 4 
relative to Ca 3 . The CCSD(T) quadratic force constants and harmonic frequencies are in 
excellent agreement with the respective MRCI quantities. Indeed, u>i and u >2 only differ 
by 1 cm -1 and w 3 differs by less than this. These comparisons suggests that the CCSD(T) 
level of theory is closely approaching the ra-particle limit for Ca 4 . 

The CCSD and CCSD(T) harmonic frequencies obtained with the larger [6s 5p 2d If] 
ANO basis set demonstrate the importance of using large one-particle basis sets in order 
to obtain highly accurate harmonic frequencies. The difference between the CCSD(T) 
harmonic frequencies in the two basis sets used is larger than the differences due to cor- 
relation treatment among the small basis results. As with Ca 3 , it is expected that the 
CCSD(T)/[6s 5 p 2d If] quadratic force constants and harmonic frequencies should be 
very reliable. 

Table 6 contains the complete cubic and quartic force field of Ca 4 obtained at the 
CCSD(T) level of theory with the larger ANO basis set. The resulting anharmonic con- 
stants are given in Table 7, while the fundamental vibrational frequencies and IR intensities 
of Ca 4 are presented in Table 8. As with Ca 3 , the absolute anharmonic corrections for 
the vibrational modes of Ca 4 are relatively small at only 3 cm 1 , 1 cm and 1 cm for 
Wl| u, 2 and u> 3 , respectively. The percentage effect on the fundamental frequencies is also 
similar to that found for Ca 3 . The IR intensity of the t-z vibration is 2.1 km/mol, which 
is substantially larger than the IR intensity of the e' vibration in Ca 3 , but again the best 
method to obtain fundamental frequencies from experiment may be an indirect approach. 


D. Summary of CCSD(T) Results for Small Be, Mg and Ca Clusters 

A summary of the equilibrium structures, vibrational frequencies, infrared intensities 
and binding energies for the alkaline-earth trimers is presented in Table 9 and for the 
tetramers in Table 10. The beryllium and magnesium cluster results are taken from our 
previous studies. 11,12 In all cases, the results are those obtained with the largest ANO 
basis set used in the particular investigation. For the trimers, Table 9 contains the MRCI 
equilibrium bond distance, harmonic frequencies and dissociation energy. The fundamental 
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frequencies were obtained by adding the CCSD(T) anharmonicity to the MRCI harmonic 
frequencies. The IR intensities were determined with the CCSD(T) method. For the 
tetramers, all of the results were obtained at the CCSD(T) level of theory since it was not 
possible to use the MRCI method in conjunction with the larger ANO basis sets for the 
tetramers. Thus the values summarized in Tables 9 and 10 represent the best computed 
quantities to date for the alkaline-earth trimers and tetramers. 

Ex ami nation of the binding energies in Table 9 indicates the expected trend based on 
bulk binding energies — that is the binding energies decrease in the order Be 3 > Ca 3 > Mg 3 . 
The vibrational frequencies, on the other hand, decrease in the order Be 3 > Mg 3 > Ca 3 , 
but this is determined in large part by mass effects, since the symmetry internal coordinate 
quadratic force constants for Mg 3 are smaller 12 than those for Ca 3 . The IR intensity of 
the e 1 mode is small for all of the alkaline-earth trimers. 

Comparison of the D e values of the tetramers shows the same trend as observed 
for the trimers, although the ratios are somewhat different. The binding energy of Be 4 
is significantly larger than that of Be 3 , though the Be 4 equilibrium bond distance is only 
0.28 ao shorter than the Be 3 value. The ratio D e (M 4 )/-D e (M 3 ) is largest for M — Mg 
and it is therefore not surprising that Mg 4 exhibits the largest reduction in bond length 
(0.50 ao) relative to the trimer. The reduction in the Ca 4 bond length (relative to the 
trimer) is about the same as that observed for Be even though the ratio Z) e (M 4 )/D e (M 3 ) 
is significantly larger for M = Be than for M = Ca. This last observation is probably due 
to the fact that the valence orbitals of Ca are larger than those of Be. 

The harmonic frequencies of the tetramers exhibit the same trend as observed for the 
trimers. Again, the Mg 4 vibrational frequencies are higher than the analogous Ca 4 values 
because of the mass effect. The fundamental frequency i/ 3 for Be 4 has a large positive 
anharmonic correction that is not observed in either Mg 4 nor Ca 4 . This Be 4 phenomenon 
was explained in some detail previously, 11 and relies on a symmetry argument applicable 
to tetrahedral X 4 species. Its apparent inapplicability to Mg 4 and Ca 4 is probably due 
to several factors, including the much weaker bonds present in Mg 4 and Ca 4 relative to 
Be 4 . The IR intensities of the tetramer t 2 vibrations are much larger than those calculated 
for the trimers. Nevertheless, the Mg 4 and Ca 4 intensities remain very small and it is 
likely that the best prospect for experimental determination of these frequencies will be an 
indirect method like negative ion photodetachment. On the other hand, the IR intensity 
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of the mode of Be 4 is certainly large enough to allow direct experimental observation 
provided that an experiment can be designed which will produce enough Be 4 . 

Conclusions 

The CCSD, CCSD(T) and MRCI electron correlation methods have been used to 
investigate the equilibrium structures, vibrational frequencies and binding energies of the 
Ca 3 and Ca 4 metallic clusters. In agreement with our earlier studies of the analogous Be 
and Mg clusters, it is found that the CCSD(T) method reproduces the MRCI r e , harmonic 
frequencies and D e values very well. The agreement between these methods is somewhat 
better for Ca 4 , but is still very good for Ca 3 . Complete cubic and quartic force fields 
of Ca 3 and Ca 4 have been determined with the CCSD(T) method in conjunction with a 
large ANO basis set and have been used to evaluate the anharmonic corrections needed to 
compute the fundamental frequencies. The anharmonic corrections have been determined 
via second-order perturbation theory. The absolute value of the anharmonic corrections 
is relatively small, although as a percentage relative to the fundamental frequencies they 
are similar to those observed previously for the Mg clusters. In spite of the fact that Ca is 
larger and more polarizable than Mg, and that Ca clusters are more strongly bound than 
Mg clusters, the IR intensities of the e' mode of Ca 3 and of the 1 2 mode of Ca 4 are small 
and similar to the analogous Mg quantities. It is unlikely that direct observation of these 
fundamentals will be possible. 

The MRCI and CCSD(T) equilibrium structures, vibrational frequencies and binding 
energies of the alkali metal (Be, Mg and Ca) trimers and tetramers have been summarized. 
The binding energies of the trimers and tetramers follow the bulk metal binding energies 
although the ratios of the small cluster D e values do not agree with the bulk metal ratios. 
The vibrational frequencies follow a different trend as the Mg n (n=3,4) frequencies are 
larger than the respective Ca values, but this is due to the larger mass of the Ca atom, 
since the symmetry internal coordinate force constants (which are independent of mass) 
for the Ca trimer and tetramer are larger than the respective Mg quantities. 
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Table 1 


Total energies (J B^), bond lengths (a 0 ), rotational 
constants (MHz) and binding energies (kcal/mol) for Ca 3 and Ca 4 . 




Basis 

E a 

r e 

A 

B 

D e 

Do 

Ca3 

CCSD 

[65 5 p 2d If] 

0.370728 

8.097 

1378 

689 

6.4 

6.1 


CCSD(T) 

[6s 5 p 2d If] 

0.378316 

7.884* 

1453 

726 

11.2 

10.8 


MRCI 

[6s 5 p 2d If) 

0.379234 

7.874 

1457 

728 

12.1 

11.7 

Ca4 

CCSD 

[5s 4 p Id] 

0.168123 

7.935 

717 


14.3 

13.6 


CCSD(T) 

[5s 4 p Id] 

0.180604 

7.806 

741 


22.2 

21.4 


MRCI 

[5s 4 p Id] 

0.179059 

7.824 

738 


22.4 

21.7 


CCSD 

[6s 5 p 2d If] 

0.180579 

7.694 

763 


20.9 

20.1 


CCSD(T) 

[6s 5 p 2d If) 

0.197496 

7.591 c 

784 


31.5 

30.6 


a The energy for Ca 3 is reported as —(£'+2030) and for Ca 4 as —(£+2707). 
b Vibrationally averaged bond lengths: r g = 7.917 ao and r a = 7.915 a$. 
c Vibrationally averaged bond lengths: r g = 7.615 ao and r a = 7.613 ao. 
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Table 2 

Symmetry internal coordinate force constants (aJ/A 2 ) 
and harmonic frequencies (cm -1 ) for Ca3. 



F n 

F 22 

w i (ai ) 

W2(e') 

CCSD 

0.04415 

0.08141 

75 

72 

CCSD(T) 

0.06952 

0.10835 

94 

83 

MRCI 

0.07126 

0.11460 

95 

85 
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Table 3 

Non-zero cubic (aJ/A 3 ) and quartic (aJ/A 4 ) force 
constants and the anharmonic constants (cm -1 ) for Ca 3 . 


-0.133 

F\ 2 a 2 a = Fi 2 b 2 b —0.152 

F2a.2a.2a = ~F2 a 2b2b —0.105 

Fun 0.156 

Fll2a2a = Fn2i26 0.163 

F\2a2a2a = — Fi2a2b26 0.116 

F 2 a 2 a 2 a 2 a = F 2 i,2fe2b26 = 3F 2a 2a2626 0.239 


Anharmonic Constants 


x li —0.60 

x 2 i -1.38 

x 2 2 —0.41 

922 0.18 
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Table 4 


Comparison of the CCSD(T) harmonic and fundamental frequencies 
of Ca 3 (cm -1 ). Infrared intensities (km/mol) are also included. 


Mode 

u 

V 

LJ — V 

I 


94 

92 

2 

0 

e' 

83 

81 

2 

0.4 
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Table 5 

Symmetry internal coordinate quadratic force constants (aJ/A 2 ) 
and harmonic frequencies (cm -1 ) for Ca 4 . 



Basis 

F n 

F22 

F 33 

^i(ai) 

<^2(e) 

^ 3 (^ 2 ) 

CCSD 

[5s 4 p Id] 

0.05350 

0.12546 

0.08688 

95 

73 

86 

CCSD(T) 

[5«s 4 p Id] 

0.07166 

0.14253 

0.10439 

110 

78 

94 

MRCI 

[5s 4 p Id] 

0.06960 

0.14091 

0.10397 

109 

77 

94 

CCSD 

[6s 5 p 2d If] 

0.07463 

0.15725 

0.11282 

113 

82 

98 

CCSD(T) 

[6s 5 p 2d 1/] 

0.09446 

0.17250 

0.13062 

127 

86 

105 
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Table 6 

Non-zero cubic (aJ/A 3 ) and quartic (aJ/A 4 ) force 
constants for Ca<j. 


*111 


- 0.110 

* 12 a 2 a — *12626 


-0.140 

*13x3x = * l$y3y = *13z3z 


-0.126 

* 2 a 2 a 2 a = — * 2 a 2626 


-0.074 

*2a3z3z — 2*2a3x3x 2*2a3y3y y/g *263x3® 

*263y3y 

-0.190 

*3z3y3z 


-0.013 

-*1111 


0.077 

* 1 12 a 2 a = *112626 


0.091 

*113x32 = *113y3y = *113z3z 


0.076 

* 12 a 2 a 2 a = * 12 a 2626 


0.058 

*1 2 a 3 z 3 z = — 2*i2a3x32 = ~ 2*12a3y3y = ^*1263z3x = 

^/|*1263y3y 

0.119 

*13 x3y3z 


- 0.001 

* 2 a 2 a 2 a 2 a = *26262626 = 3*2a2a2626 


0.162 

*2a2a3z3z 


0.176 

*26263z3z 


0.044 

a * 2 a 2 a 3 y 3 y = * 2 a 2 a 3 x 3 x = |(*2a2a3z3z + 3*26263z3z) 


0.077 

a *26263y3y ^ *26263x3x = 4 (3*2a2a3z3z + *26263z3z) 


0.143 

*2a263y3y — *2a263x3x 4 ( *2 a 2 a 3 z 3 z *26263z3z) 


0.057 

*323x3x32 — *3y3y3y3y *3z3z3z3z 


0.262 

*3x3x3y3y = *3x3x3z3z — *3y3y3z3z 


0.043 


“ Dependent force constants related to F2a.2a.zzZz and F2b2bZzZz- 
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Table 7 

Anharmonic constants (cm -1 ) for Ca 4 . 


Anharmonic Constants 


111 

-0.31 

Z 21 

-0.52 

I 22 

-0.13 

2:31 

-0.83 

^32 

-0.35 

Z 33 

-0.18 

9 22 

0.09 

9 33 

0.03 

1 23 

-0.05 

^33 

-0.03 
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Table 8 

Comparison of the CCSD(T) harmonic and fundamental frequencies 
of Ca 4 (cm -1 ). Infrared intensities (km/mol) are also included. 


Mode 

C 0 

V 

U) — V 

I 

aa 

127 

124 

3 

0 

e 

86 

85 

1 

0 

^2 

105 

104 

1 

2.1 
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Table 9 

Su mm ary of results for the alkaline earth trimers. 0. 



r e 

wi(a'i) 

«a(e') 

V\ 

^2 

De 

Be 3 

4.200 

490 

427(0.5) 

469 

410 

22.5 

Mg 3 

6.373 

110 

115(0.2) 

101 

109 

6.3 

Ca 3 

7.874 

95 

85(0.4) 

93 

83 

12.1 


® Units are a 0 for r e , cm -1 for the harmonic and fundamental frequencies, and kcal/mol 
for D e . The value in parentheses is the IR intensity in km/mol. The Be 3 results are from 
reference 11 and the Mg 3 results are from reference 12. See text for details of correlation 
methods used. 
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Table 10 

Su mm ary of CCSD(T) results for the alkaline earth tetramers. a 



r e 

uq(ai) 

w 2 (e) 

013 (^ 2 ) 

Vi 

^2 

^3 

D e 

Be 4 

3.921 

663 

469 

571(29.7) 

639 

455 

682 

79.5 

Mg 4 

5.877 

192 

147 

171(2.4) 

184 

143 

167 

23.9 

Ca 4 

7.591 

127 

86 

105(2.1) 

124 

85 

104 

31.5 


a Units are a 0 for r e , cm -1 for the harmonic and fundamental frequencies, and kcal/mol 
for D e . The value in parentheses is the IR intensity in km/mol. The Be 4 results are from 
reference 11 and the Mg 4 results are from reference 12. 
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Abstract 

We show that by a suitable change of variables, the derivatives of molecular 
integrals over Gaussian-type functions required for analytic energy derivatives can 
be evaluated with significantly less computational effort than current formulations. 
The reduction in effort increases with the order of differentiation. 


I. Introduction 

Analytic energy derivative methods have revolutionized the application of com- 
putational quantum chemistry to problems of chemical interest [1]. The location 
and characterization of stationary points on polyatomic molecular potential energy 
surfaces can be accomplished so much more efficiently using analytic derivatives 
than with techniques based on computing energies alone that the development and 
extension of analytic derivative methods has been one of the most active fields of 
methodological research in quantum chemistry in recent years. Given the gradi- 
ent and Hessian of the energy with respect to the nuclear coordinates, a variety of 
strategies have been developed that are guaranteed to converge to minima on po- 
tential surfaces and that can efficiently locate other stationary points, particularly 

t Mailing address: NASA Ames Research Center, Moffett Field, California 94035- 
1000 
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transition states. These strategies can also be used to “walk” on surfaces from one 
minimum to another, thereby defining a reaction coordinate, and among the most 
elegant and conceptually illuminating studies of this sort are the investigations 
of Ruedenberg and co-workers on rearrangement reactions of small hydrocarbon 
species (see Refs. 2-5 and references therein). It is thus a great pleasure to dedicate 
this contribution to Professor Ruedenberg on the occasion of his 70th birthday. 

Of course, in order to perform such walks and optimizations it is imperative 
to evaluate the energy derivatives efficiently at the computational level of interest 
(Hartree-Fock or some correlated treatment). As noted above, much work has been 
performed in this area, and several reviews are available [1,6,7]. We shall concentrate 
here on a topic that ultimately affects the computational effort necessary to evaluate 
energy derivatives for any ab initio method that relies on a basis set expansion of 
Gaussian one-electron functions. 

Wave functions for polyatomic molecules are invariably expanded in a basis set 
that is centred on the various nuclei, and so in a calculation of the energy derivative 
of nth order with respect to the nuclear coordinates, up to nth-order derivatives of 
the one- and two-electron integrals are required. These derivative integrals can in- 
volve differentiation of the operators as well as differentiation of the basis functions, 
but the greatest computational problems arise from the differentiation of the basis 
functions. Like the evaluation of integrals over Gaussians [8,9], the calculation of 
integrals over differentiated Gaussians has been the subject of many investigations 
and numerous efficient computational schemes have been devised. In this work we 
show how the efficiency of derivative integral evaluation can be improved by some 
simple manipulations of variables. We shall briefly review the McMurchie-Davidson 
scheme [8] for computing Gaussian integrals and derivative integrals, and then show 
how a change of differentiation variables simplifies the formulas. 

II. Derivative Integral Formulas 

We shall expand the Gaussian charge distributions that appear in the integrals 
in Hermite functions, as described by McMurchie and Davidson [8] (see also Saun- 
ders [9]). Let us represent an unnormalized Cartesian Gaussian function centred 
at A by 

G ijk {r, a, A) = x^y^z'X exp(-ar^), (1) 

where x A = x- A x , etc. We can consider one Cartesian direction, say x, represented 
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as 

Gi{x,a,A x ) = x l A exv{-ax 2 A ). 

The overlap distribution of two such functions is expanded as 

a, 6 , A x , B x ) = Gi(x,a,A x )Gj(x,b,B x ) 

*+ j 

= Y^E\ 3 {a,b,A x ,B x ) A t (x,p,P x ), 
t= o 

where the Hermite function A.t(x,p, P x ) is defined by 

A t (s,p,Px) = (d/dP x y exp (-px 2 p) 

with , 

P P 

and 

p = a + 6. 

The expansion coefficients 6, u4 x , are obtained from 

K 

where 


i+1,i = - -R x El j + (f + 1 )£*% 

2p P 


and 


— A x B x 

Eq° = exp (-—Rl). 

V 


( 2 ) 


( 3 ) 

( 4 ) 

(5) 

( 6 ) 

( 7 ) 

( 8 ) 
( 9 ) 


Henceforth we shall not always fist the arguments of the expansion coefficients or 
Hermite functions, but we wish to emphasize here that the expansion coefficients 
depend on u, 6, and R x only, while the Hermite functions are independent of R x . 


i+j 


Sl ij {z,a,b,A x ,B x ) = J^ E * i ( a , b , R t ) A t (x,p,P x ). 


( 10 ) 




In terms of the Hermite functions and expansion coefficients we can express a 
two-electron integral 


// 


*aVa*a ex P(~ ar A ) x bVb z b exp(-6r|) 

x r ~2 XcVc z c‘ e x v(- cr c) x dV 1 D z D exp (-dr 2 D )dr l dr 2 (11) 
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as 


H-j 


i'+3 


Y,El j {a,b,A x ,B x ) Y E V ( c ,d,C x ,D x ) 

t=0 t' = 0 

k+l k'+l‘ 

xY E “( a ’ b ’Ay,By) Y E u' l \ c > d > c v> D v) 


U-O 
m-\- n 


u* = 0 


x y K n { a i b ’ A *i B *) Y E ™‘ ( c , d ’ c *’ D *) 

t )=0 

x (iuvlr^ 1 1 t'u'v 1 ), 


v'—O 


( 12 ) 


where 

(tuv\r^2 | t*u v) 

A t (x,p, P x ) At'(z, g, Q x ) Au(2/j?5 Py) Au'(y5 ?5 Qy ) 

X A v {z,p,P z )Av‘(z,q,Qz)r^2 drjdr 2 . ( 13 ) 

and q and Q are defined analogously to p and P but for the second charge dis- 
tribution. Thus in practice we evaluate integrals over the Hermite function basis 
and combine those with the expansion coefficients to give integrals over primitive 
Gaussians. Some modifications to the form (12) are desirable from the point of 
view of efficiency, as discussed by Saunders [9], but for schematic purposes we can 
use (12). The first step, evaluation of the Hermite function integrals, is fast. The 
second step, which we can regard as a transformation from the Hermite function 
basis to the Cartesian Gaussian basis, is relatively time-consuming and is certainly 
more expensive than calculating the Hermite function integrals. Finally, if required, 
we combine these integrals with basis set contraction coefficients to give the final 
integrals. In fact, some of the expansion steps can be taken outside the contraction 
step, with a consequent improvement in efficiency. 

In a derivative integral we are interested in derivatives of Cliji dflij / dA x and 
dClij/dB x for first derivatives, for example. Conventionally, we would differentiate 
the orbitals (2) first and then expand the overlap distributions of the differentiated 
orbitals analogously to fly above. For example, for the derivative with respect 
to A x we obtain 

dClij 


i+j+l 

Y F t A <- (14) 
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Note that the sum here is over more terms than appear in the undifferentiated charge 
distribution (3) — higher orders of differentiation would increase this summation 
range further. The new coefficients are defined in terms of the coefficients E t 
above by 

Fi j = 2aE l t +i ' j - iE\~ l ' j . (15) 


Analogous coefficients can be defined for higher orders of differentiation or for dif- 
ferentiation with respect to B x . In this approach, then, we compute derivative inte- 
grals using the same general scheme (12) as for undifferentiated integrals. Since the 
expansion of the differentiated charge distributions in Hermite functions is longer 
than for the undifferentiated distributions, the work required to transform from the 


Hermite function basis to the Cartesian Gaussian basis is greater. Further, as the 
order of differentiation increases this extra work becomes larger and larger. Hence 
this approach is not well-suited to higher derivatives. 

Let us instead consider differentiation with respect to the variables P x and R Xi 

for which 


and 


d 

a d 

d 

“f" ^ 

(16) 

dA x ' 

p dP x 

dR x 

d 

b d 

d 

(17) 

dB x ‘ 

V dP x 

dR x ’ 


We recall that the Hermite functions are independent of R x , while the expansion 
coefficients are independent of P x . Hence we can expect the expressions for the 
differentiated charge distributions to be simpler in terms of these variables, although 
we must eventually transform the derivatives back to the A X ,B X representation. We 
obtain for the derivatives 


dClij _ yi „ij dA. t 

~ h ' dp ■ 


£ E? A, +1 

i=0 


(18) 


CbJ-LVJ. . , . . , 

dClij _ dE^ 1 

™*~k dR * 

Denoting dE\ j jdR x by E ^' 1 , we obtain the expansion relation 

= ±-E*\ - -(RxE 1 /’ 1 + E?) + (t + 1)E%\ 

Zp p 


(19) 


( 20 ) 
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by differentiating the relation (7) above. 

We can make several important observations about these derivative formulas. 
First, the combination of expansion coefficients and Hermite functions in (18) above 
is over exactly the same range as the summation to give undifferentiated integrals: 
the only difference is that the degree of the Hermite function has increased by 
one. Hence the code required to evaluate this term is the same as required in the 
undifferentiated case, and the number of operations is also the same. (It is easy 
to see that this holds true in any order of differentiation for this term.) As we 
saw above, this is not the case if we differentiate with respect to the variables A x 
and B x , because then a linear combination of different degree Hermite functions 
and expansion coefficients appears. 

Second, calculation of the differentiated expansion coefficients E t ' requires 
essentially the same code again as for the undifferentiated case, with the obvious 
addition of an extra term in the expansion relation, and a starting value 

( 21 ) 

V 

obtained by differentiating (9). As we noted, the index range of the coefficients that 
are required is the same as that for the undifferentiated case, so the actual work 
required to combine Hermite function integrals and expansion coefficients does not 
increase. (The precomputation of the expansion coefficients themselves is of course 
a very rapid step.) 

Third, in the usual scheme the index range of the program loops over the 
variables t, u, v depends on the direction of differentiation (i.e., differentiation with 
respect to A x , A y , etc). Thus these loops must be executed with different ranges 
for each of the three directions for first derivative integrals, for example. With 
our transformation of variables, the loop index ranges become independent of the 
direction of differentiation, so the program logic is simplified and the overheads 
are reduced. We may also note here that this approach in no way diminishes the 
possibilities for vectorizing the calculation of the integral derivatives. Indeed, the 
simplifications to the program loop structure are likely to enhance these possibilities. 

Fourth, we can obtain an additional simplification as follows. Adding (16) 
and (17) we obtain 

_JL = A d — ( 22 ) 

dB x dP x dA x ' 
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Now, (in addition to saving one multiplication) this form of the expression for 
the derivative with respect to B x does not depend on the orbital exponents at 
all. Hence we can delay the transformation to the B x derivative until later in the 
calculation, for example, until after the contraction step, so that the time required 
for this variable transformation becomes negligible. This is most important for first 
derivatives, as in any order of differentiation only one term can be treated this way. 

In the case of higher derivatives there is a variety of terms to be considered but 
the scheme remains essentially the same. For example, the nth-order differentiated 
expansion coefficients with respect to R x are obtained from the recursion formula 

Et +1 ' j ' n = — E*? - -{R*E? n + nEf' 1 ) + (t + l)Ei$ (23) 

2 p P 

with starting values 

E?» in+I = -—(} T” -1 ) < 24 ) 

V 

and the identification 

E l t j '° = E l t j ( 25 ) 

Higher derivatives of the Hermite functions with respect to P x (19) are trivially 
obtained. We note further that if the two charge distributions that appear in an 
integral are differentiated separately, the total savings is the product of the individ- 
ual reductions in work, since the two differentiations are independent. For multiple 
differentiation of the same charge distribution, we recall that by using our trans- 
formation of differentiation variables the summation range in the Hermite function 
to Cartesian Gaussian transformation is independent of the order of differentiation. 
Hence the savings increase as the order of differentiation increases, since in the 
conventional scheme the work required to accomplish this transformation increases 
substantially with the order of differentiation. In order to obtain an estimate of 
what savings are possible, we must also include an estimate of the effort required 
to transform back to the A X ,B X representation. We shall now present operation 
counts showing that it is always preferable to use our transformation of differenti- 
ation variables. 

In order to simplify the counting we consider only floating-point operations 
(multiplication and addition), which are weighted equally. In addition, in our count 
we have not taken advantage of the possibility of deferring transformation of some 
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derivatives until after contraction: in effect, we are counting operations only for 
primitive Gaussians and ignoring any additional savings that might accrue from 
moving manipulations outside the contraction step. If anything, neglecting this 
possibility favours the conventional approach to derivative integrals. 

We have listed operation counts for differentiation of SS , PP , and DD distri- 
butions in Table 1. We have not included the calculation of the Hermite function 
integrals, which is fast and contributes the same work to both cases, the conventional 
approach and our new scheme. Further, the transformation of the second charge 
distribution in the integral has also been excluded. We see that for the SS case the 
total operation count is not much affected by whether or not the transformation of 
variables is performed. However, for higher angular momentum functions there is 
a decided advantage to using the transformation of variables, and this advantage is 
clearly growing with the order of differentiation. As a further illustration of this, we 
note that for third derivatives of a PP distribution, for example, the conventional 
method would require 14 448 operations, while using the transformation of variables 
the work would be reduced to 8 340 operations: a savings of 42%. 

Finally, some other aspects of this scheme deserve comment. We note that 

9 = b d _ a d (26) 

dR x p dA x p dB x 


Therefore, the operation is not the same as the differentiation QA ^ — But 

if A and B coincide then the differentiation with respect to R x does not contribute 
to the energy derivative: only the differentiation with respect to P x contributes. 
This simplification is already used in the ABACUS program [10]. We also note 
that the use of translational invariance to reduce the computational labour is not 
affected by our transformation of variables: for first derivatives, for example, we 
have 

(27) 


dl dl 

1 = 0 , 

dP x dQ x 


where I represents the two-electron integral in (11), from the use of translational 
invariance. 


Conclusions 

We have shown that by employing a transformation of differentiation variables, 
the work required to evaluate derivative integrals can be substantially reduced. The 
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advantages of our new approach increase both with the order of differentiation and 
with the angular momentum of the Gaussian functions involved. Savings will be 
obtained in the calculation of energy derivatives for any wave function that is ex- 
panded in a Gaussian basis. In particular, the economies obtained by applying these 
methods to the calculation of third or higher derivative integrals will be substantial. 
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Table 1. Operation counts for differentiation. 



ss 

PP 

DD 

First 

Derivatives 



Hermite/ Cartesian Transformation 

12 

396 

4 032 

P X ,R X to A X ,B X Transformation 

9 

81 

324 

Total 

21 

477 

4 356 

Conventional 

24 

672 

6 144 


Second Derivatives 


Hermite/ Cartesian Transformation 

42 

1 386 

14 112 

P X ,R X to A X ,B X Transformation 

93 

837 

3 328 

Total 

135 

2 223 

17 460 

Conventional 

150 

3 678 

30 912 
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Abstract 


We have calculated the second hyperpolarizability 7 of neon using the CCSD(T) 
method. The accuracy of the CCSD(T) approach has been established by explicit compar- 
ison with the single, double and triple excitation coupled-cluster (CCSDT) method using 
extended basis sets that are known to be adequate for the description of 7 . Our best es- 
timate for 70 of 110±3 a.u. is in good agreement with other recent theoretical values and 
with Shelton’s recent experimental estimate of 108±2 a.u. Comparison of the MP2 and 
CCSD(T) hyperpolarizability values indicates that MP2 gives a very good description of 
the electron correlation contribution to 70 . We have combined MP2 frequency-dependent 
corrections with the CCSD(T) 70 to yield values of 7 (- 2 u;;u;,a>, 0) and 7 K (-w;u;, 0 , 0). 


2 



1. Introduction 


Theoretical determination of hyperpolarizabilities has been a topic of much interest 
recently, since knowledge of atomic and molecular hyperpolarizabilities is central to the 
understanding of the non-linear response of matter to light. In particular, organic materials 
with large hyperpolarizabilities are candidates for applications such as optical switching 
and second harmonic generation, and there is great potential for interaction between theory 
and experiment in the study of these systems. 

From a theoretical point of view, it is important to understand the requirements for 
determining accurate hyperpolarizabilities for small systems, because it is possible to use 
large one-particle basis sets and sophisticated electron correlation treatments for these 
species, and thereby to evaluate the effects of approximations that will be necessary for 
the study of the hyperpolarizabilities of larger systems. Hence for small systems it is 
desirable to estimate the accuracy of the calculated hyperpolarizability. This may be 
accomplished in two ways. First, the quality of the one and n-particle approximations 
used in the calculation can be systematically improved and the convergence of the result 
can be monitored. This is perhaps the preferred approach from a theoretical standpoint. 
Alternatively, the theoretical value can be compared directly with experiment, although 
the possibility of error cancellation between the one and n-particle approximations must 
always be borne in mind. 

Study of the hyperpolarizabilities of the rare gas atoms has a number of advantages. 
In particular, for neon sophisticated levels of theory and large one-particle basis sets can 
be employed. Experimental gas-phase electric-field-induced second harmonic generation 
data for the rare gases is available over a range of frequencies 1 , which makes extrapolation 
to the static limit possible for the purposes of comparison with a theoretical static value. 
In addition, vibrational effects (which have been shown to be non-negligible for some 
molecular values, see, for example, Refs. 2 and 3) vanish for atoms. 
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It must be noted, however, that it is a non-trivial task to demonstrate convergence 
of calculated values for hyperpolarizabilities. Previous work has shown that there can be 
a much larger electron correlation contribution to hyperpolarizabilities than for linear po- 
larizabilities, and that the contribution of higher excitations is not insignificant (see, for 
example, Refs. 4-6). For example, in the case of 7 of neon, where the effect of electron 
correlation is about 40 a.u., or 40% of 7, the perturbational estimate of triple excita- 
tions contributes 8 a.u. or 20% of the total electron correlation contribution. 4 Since the 
contribution of higher excitations is so large, the applicability of approximate methods 
for estimating the effects of higher excitations in hyperpolarizability calculations might 
be questioned. For example, the single and double excitation coupled-cluster method in- 
cluding an estimate of triple excitations through the fourth and partially the fifth order 
of perturbation theory, (CCSD(T)) 7 has had great success in describing the structure and 
frequencies of a number of ‘difficult’ chemical systems — that is, systems whose wave 
functions are not dominated by a single reference. 8 ' 9 This success notwithstanding, it is 
essential to investigate the utility of this approach specifically for determining accurate 
hyperpolarizabilities. 

The reliability of a correlation treatment is best evaluated by comparison with a full 
configuration-interaction (Cl) calculation in the same one-particle basis set. However, since 
even at the self-consistent field (SCF) level of theory diffuse / type functions contribute 
10 a.u. to the hyperpolarizability of Ne 4 , a full Cl calibration in a realistic basis set is not 
feasible. Here, in order to establish the accuracy of the CCSD(T) method for 7 of neon, 
we compare instead to results obtained with the full single, double and triple excitation 
coupled-cluster method (CCSDT), 10 ' 11 in a basis set which is known to be adequate for 
the description of hyperpolarizabilities. We can thus assess the accuracy of the computed 
CCSD(T) value, as well as compare it with experimental and other theoretical values. We 
note that a previous study has demonstrated that for correlation energies, CCSD(T) is an 
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excellent approximation to CCSDT. 12 


2. Computational Methods 

The one-particle basis sets used in this work are similar to those employed previously. 4 
They were derived from van Duijneveldt’s (13s8p) primitive set 13 augmented with a (6d4/) 
polarization set, with exponents chosen as an even-tempered sequence <*= 2.5 n a 0 -, n 
— o,...,fc with a 0 = 0.20, 0.61 for the d and / functions, respectively. This was con- 
tracted to [4s 3p 2d If ] using atomic natural orbitals. 14 The two outermost sets of 
spd functions and the outermost / function were uncontracted to give basis C, denoted 
[ 444 - 1 - 1,8 3 _| 4 -pip 2+1+ld 1+1/ ] • We use the notation C in order to be consistent with 
our previous study. 4 Additional diffuse functions were then added by extrapolating from 
the outermost function in an even-tempered sequence, <*= 2.5 _n oio- For example, the ad- 
dition of one set of diffuse functions is denoted + (lslpldlf). In some calculations the 
basis was further augmented with two diffuse g functions with exponents of 0.29 and 0.11. 
Basis set C was completely uncontracted, and two tight d functions were included (a«i 
= 123.53, 49.41) for calculations in which core correlation was included. Only the true 
spherical harmonic components of the basis functions were used throughout. 

Energies were calculated using self-consistent field (SCF), single and double excitation 
coupled-cluster (CCSD) and second-order Mpller-Plesset perturbation theory (MP2). The 
effect of triple excitations was investigated using (a) the CCSD(T) method, which includes 
an estimate of the triples through fourth and partially fifth order of perturbation theory , 
based on the CCSD amplitudes in the perturbation energy expressions; (b) the CCSDT 
method, which explicitly includes all single, double and triple excitations; and (c) the 
CCSD + T(CCSD) method 15 which includes only the fourth-order perturbation theory 
contribution based on CCSD amplitudes. 

The dipole polarizabilities are defined 16 by the energy response to an applied electric 
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field of strength F: 


E(F) = Eo - \aF> - 

Electric fields of 0, 0.002, 0.004, 0.008, 0.012, 0.016 and 0.020 a.u. were applied and 
the energy responses were fitted to the sixth-order polynomial in the field strength. The 
SCF value of 7 obtained from the fit agrees with the SCF value obtained from finite 
displacements of analytic j3 values to within 0.1 a.u. Tests of the fit for the correlated 
values indicate the error in 7 due to the fit is less than 0.2 a.u. 

The SCF, MP2 and most of the coupled-cluster calculations were performed using 
the MOLECULE-SWEDEN , 17 CADPAC, 18 , and TITAN 19 programs. The CCSDT calcu- 
lations were performed with a program written by one of us (GES ). 11 The SCF energies 
were converged to 10 18 E ^ or better and the CCSD and CCSDT energies to 10 

3. Results and Discussion 

The values for the linear polarizability, a and the hyperpolarizability, 7 determined 
in this work are summarized in Table 1. We note first that the SCF, MP2, CCSD and 
CCSD(T) results for both a and 7 are essentially identical in basis sets C+(3s3p2d3/) and 
C+(2s2pld2/). This establishes that these values are converged with respect to further 
addition of diffuse s, p, d,and / functions in the one-particle basis set. It also verifies 
that basis C+(2s2pld2f) is a good choice for comparison of the CCSDT and CCSD(T) 
results. The CCSDT 7 value of 110.9 a.u. in this basis is in very good agreement with 
that from the CCSD(T) method ( 111.2 a.u.). This establishes that the (T) correction is 
a reliable estimate of the contribution from connected triple excitations to the hyperpo- 
larizability of Ne, even though the contribution is large. Comparison of the CCSD(T) 
and CCSD + T(CCSD) results with the C+(Zs3p2d3f) basis shows that the fourth-order 
correction alone is a substantial overestimate. 

The effects of diffuse higher-order angular momentum functions and of core correla- 
tion have also been investigated at the CCSD(T) level of theory, in order to improve upon 
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our previous best computed value of 111 a.u . 4 and to assess further our error estimate of 
±4 a.u. Uncontraction of the one-particle basis set reduces 7 from 111.0 a.u. to 110.3 a.u., 
and core correlation gives a further reduction of 0.5 a.u. Conversely, diffuse g functions 
increase 7 by 0.2 a.u. These changes are very similar to the values of —0.7, —0.5, and 
0.2 a.u., respectively, obtained at the CCSD level of theory, whereas the MP2 values indi- 
cate that second-order perturbation theory may underestimate the effects of uncontraction 
(_ 0.4 a.u.) and core correlation (—0.3 a.u.). This is consistent with the results for argon, 
where the CCSD estimate of the reduction in 7 due to core correlation is larger than the 
MP2 value. The CCSD(T) results thus support our earlier conclusions 4 that the combina- 
tion of further augmentation of the one-particle basis set and further improvement in the 
treatment of core correlation is likely to have little effect on 7 * 

The importance of higher than triple excitations also requires discussion. Single and 
double excitations increase 7 by 34 a.u. and triple excitations by a further 9 a.u. However, 
the coupled-cluster expansion should be rapidly convergent, since all the corresponding 
disconnected terms are included to infinite order at each level of CC theory. The next 
most important contribution comes from connected quadruple excitations: these can be 
expected to contribute substantially less than connected triple excitations. It is highly 
unlikely that these terms contribute as much as 3 a.u. We have previously established 4 
that the CCSD(T) energy, at least, agrees almost perfectly with full Cl benchmark results 
for Ne. We thus arrive at a best estimate of 110 a.u., with an uncertainty of 3 a.u. The 
latter is almost entirely due to uncertainty in the estimated contribution of connected 
quadruple (and higher) excitations. We can compare our result with Shelton’s most recent 
estimate of 108±2 a.u . 21 based on extrapolation of electric-field-induced second harmonic 
generation values. This differs from his earlier static value of 119±2 due to his discovery 
of a systematic error in some of the measurements . 21 Our value is also in good agreement 
with other recent theoretical values. These include Maroulis and Thakkar’s estimate of 
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114±9 a.u. based on the double-excitation coupled-cluster method augmented with a 
perturbational estimate of single and triple excitations, denoted CCD+ST(CCD), 22 as 
wed as Chong and Langhoff’s CCSD(T) value of 111.0 a.u. 23 , and 108±5 a.u. from the 
restricted active space self-consistent field calculations of Jensen and co-workers. 24 

It is clearly much easier to demonstrate convergence for a than for 7. The electron 
correlation contribution is smaller (around 11%) for a and the CCSD(T) and CCSDT 
values agree to 3 decimal places. The diffuse function requirement is also less stringent 
than for 7. The largest remaining corrections are for core correlation and uncontraction 
of the one-particle basis set which decrease a by 0.013 a.u., or 0.5%. Our best computed 
value for a is 2.677 a.u. to which we conservatively assign an uncertainty of 0.015 a.u. 
This is in excellent agreement with the value of 2.669 a.u. derived from dipole oscillator 
strength distributions by Kumar and Meath. 25 

Finally, the rather good agreement between the MP2 and CC polarizability and hy- 
perpolarizability suggests that the former method may be useful even where, as is the case 
for neon, the correlation contributions are large. Because it is much simpler to evaluate 
dynamic polarizabilities with MP2 than with CC methods, this suggests that an efficient 
route to reliable frequency-dependent results may be to correct CC static values with MP2 
differences between static and dynamic values. It is certainly the case that the MP2 values 
in themselves would be much more reliable than values obtained using SCF methods to 
describe the frequency dependence. We may also note that the experimental frequencies 
are far from resonance, so an error in the frequency-dependent hyperpolarizability arising 
from the error in the prediction of the poles at the MP2 level of theory is not a matter for 

concern. 

In this work we have combined the MP2 frequency-dependent corrections 26 for 
7 (— 2w;u>,u>,0) of neon with our best estimate for the static value and arrive at 
112±3 a.u.(A=1319 nm), 113±3 a.u.(A=1064 nm), 124±4 a.u.(A=514.5 nm), to 
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be compared with experimental second harmonic generation values of 111.1±0.8 a.u. , 
109.9±0.5 a.u. 1 and 122.2±0.5 a.u. 21 , respectively. The theoretical results show no nega- 
tive dispersion (note that the uncertainties in the theoretical numbers are not independent 
and are likely to be strongly correlated with one another), and are also in good agreement 
with the revised experimental values. We may also combine the MP2 frequency-dependent 
correction 26 for 7*(-w;w,0,0) of neon with our best estimate of the static value and ob- 
tain 113±3 a.u. for A =632.8 nm. This value is somewhat higher than the Kerr effect 
measurement of 101 ±8 a.u. 27 
4. Conclusions 

We have demonstrated by comparison with full CCSDT results that the CCSD(T) 
method provides an accurate description of the hyperpolarizability of neon. Since the 
CCSD(T) method has the advantage of being rather inexpensive, it allows extensive inves- 
tigation of the one-particle basis set requirements and the core correlation treatment for 
a and 7 of neon. Our best computed CCSD(T) value for 7 of neon, including the effects 
of multiple sets of diffuse a, p, d and / functions and of core correlation, is 109.8 a.u. 
After incorporating corrections for diffuse g functions and a more complete treatment of 
triple excitations our best estimate is 110 a.u. with an error estimate of 3 a.u. This re- 
sult is in agreement with all of the most recent theoretical values 22-24 and in line with 

1 21 

the static value extrapolated from experimental frequency-dependent measurements. 
Together with our assignment of uncertainty, this fulfils our criteria for establishing the 
accuracy of the computed value. 


9 



References 


1. D.P. Shelton, Phys. Rev. A 42, 2578 (1990). 

2. D.M. Bishop, Rev. Mod. Phys. 62, 343 (1990). 

3. D.M. Bishop and B. Kirtman, J. Chem. Phys. 95, 2646 (1991). 

4. P.R. Taylor, T.J. Lee, J.E. Rice and J. Almlof, Chem. Phys. Lett. 163, 359 (1989); 
Chem. Phys. Lett, xxx, xxx (1991). 

5. R.J. Bartlett and G.D. Purvis, Phys. Rev. A20, 1313 (1979). 

6. G. Maroulis and A.J. Thakkar, J. Chem. Phys. 93, 652 (1990). 

7. K. Raghavachari, G.W. Trucks, J.A. Pople and M. Head-Gordon, Chem. Phys. Lett. 
157, 479 (1989). 

8. T.J. Lee and G.E. Scuseria, J. Chem. Phys. 93, 489 (1990). 

9. T.J. Lee, A.P. Rendell and P.R. Taylor, J. Chem. Phys. 92, 489 (1990). 

10. J. Noga and R. J. Bartlett, J. Chem. Phys. 85, 2779 (1986). 

11. G. E. Scuseria and H. F. Schaefer, Chem. Phys. Lett. 152, 382 (1988). 

12. G. E. Scuseria and T. J. Lee, J. Chem. Phys. 93, 5851 (1990). 

13. F. B. van Duijneveldt, IBM Publication RJ945 (1971). 

14. J. Almlof and P.R. Taylor, J. Chem. Phys. 86, 4070 (1987). 

15. M. Urban, J. Noga, S.J. Cole and R.J. Bartlett, J. Chem. Phys. 83, 4041 (1985). 

16. A.D. Buckingham, Adv. Chem. Phys., 12, 107 (1967). 

17. MOLECULE-SWEDEN is an electronic structure program system written by J. 
Almlof, C. W. Bauschlicher, M. R. A. Blomberg, D. P. Chong, A. Heiberg, S. R. 
Langhoff, P.-A. Malmqvist, A. P. Rendell, B. 0. Roos, P. E. M. Siegbahn, and P. R. 
Taylor. 

18. R. D. Amos and J. E. Rice, CADPAC:Cambridge Analytic Derivatives Package, Issue 
4.0 (Cambridge University, Cambridge, England, 1987). 


10 



19. TITAN is a set of electronic structure programs written by T.J. Lee, A.P. Rendell and 
J.E. Rice. 

20. J.E. Rice, P.R. Taylor, T.J. Lee and J. Almlof, J. Chem. Phys. 94, 4972 (1991). 

21. D.P. Shelton, personal communication and work presented at the North American 
Chemical Congress (ACS), August, 1991. 

22. G. Maroulis and A.J. Thakkar, Chem. Phys. Lett. 156, 87 (1989). 

23. D.P. Chong and S.R. Langhoff, J. Chem. Phys. 93, 570 (1990). 

24. H.J. Aa. Jensen, P. Jprgensen, H. Hettema and J. Olsen, Chem. Phys. Lett., in press. 

25. A. Kumar and W.J. Meath, Can. J. Chem., 63, 1616 (1985). 

26. J.E. Rice, J. Chem. Phys., submitted. 

27. A. D. Buckingham and D. A. Dunmur, Trans. Faraday Soc., 64, 1776 (1968). 


11 



Table 1 

Neon dipole polarizabilities (a.u.) 


Basis set 

Method 

a 

7 

C a + (2s2pld2 f) 

SCF 

2.377 

68.66 


MP2 

2.713 

110.6 


CCSD 

2.643 

102.2 


CCSD(T) 

2.690 

110.9 


CCSDT 

2.690 

111.2 

C a + (3s3p2d3/) 

SCF 

2.377 

68.68 


MP2 

2.713 

110.7 


CCSD 

2.643 

102.2 


CCSD + T(CCSD) 

2.703 

115.3 


CCSD(T) 

2.690 

111.0 

C a + (ZsZp2dZf2g) 

SCF 

2.377 

68.67 


MP2 

2.716 

110.9 


CCSD 

2.645 

102.4 


CCSD(T) 

2.692 

111.2 

{lZs8p8d4f)+{ZsZp2dZf) b 

SCF 

2.377 

68.67 


MP2 

2.708 

110.3 


CCSD 

2.636 

101.5 


CCSD(T) 

2.684 

110.3 

(13s8p8d4/)+(3s3p2d3/) c 

MP2 

2.703 

110.0 


CCSD 

2.628 

100.8 


CCSD(T) 

2.677 

109.8 


a See text for the definition of C. 

6 Uncontracted basis with 8 electrons correlated 
c Uncontracted basis with 10 electrons correlated 
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ABSTRACT 

We have investigated the termolecular reaction involving concerted hydrogen exchange 
between three HF molecules, with particular emphasis on the effects of correlation at the 
various stationary points along the reaction. Using an extended basis, we have located 
the geometries of the stable hydrogen-bonded trimer, which is of C 3 a symmetry, and the 
transition state for hydrogen exchange, which is of D 3h symmetry. The energetics of the 
exchange reaction were then evaluated at the correlated level, using a large atomic natural 
orbital basis and correlating all valence electrons. Several correlation treatments were used, 
namely, configuration interaction with single and double excitations, coupled-pair func- 
tional, and coupled-cluster methods. We are thus able to measure the effect of accounting 
for size-extensivity. Zero-point corrections to the correlated level energetics were deter- 
mined using analytic second derivative techniques at the SCF level. Our best calculations, 
which include the effects of connected triple excitations in the coupled-cluster procedure, 
indicate that the trimer is bound by 9±1 kcal/mol relative to three separated monomers, in 
excellent agreement with previous estimates. The barrier to concerted hydrogen exchange 
is 15 kcal/mol above the trimer, or only 4.7 kcal/mol above three separated monomers. 
Thus the barrier to hydrogen exchange between HF molecules via this termolecular process 

is very low. 


f Mailing address: NASA Ames Research Center, Moffett Field, CA 94035-1000 
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The HF trimer is known to adopt a C 3 h equilibrium structure 4 , while the transition 
state for concerted hydrogen exchange has been found to display D-$h symmetry. We shall 
consider the following three reactions: 


3HF — *• (HF) 3 [C 3h ] 

(Rl) 

(HF) 3 [C 3ft ] (HF) 3 [D 3h ] 

(it 2) 

3HF (HF) 3 {D 3h }. 

(it 3) 


Clearly, the energies of these three reactions are not independent, but it is convenient to 
retain all three for purposes of discussion. 

There have been a number of recent calculations of the relative energies of the C 3h and 
D 3 h forms of the hydrogen fluoride trimer, that is, the energy of (f?2), generally at the self- 
consistent field (SCF) level. Gaw and coworkers 4 calculated this energy to be 29.5 kcal/mol 
using a double-zeta plus polarization (DZP) basis. From their data, the actual energy for 
reaction (i?3) to proceed is 14.6 kcal/mol, uncorrected for vibrational effects. Heidrich 0 
et al. employed a variety of basis sets with geometries obtained by optimization at the 
4-31G split-valence level. Their best calculation gave AE(D 3 k-C 3 h)= 37.1 kcal/mol (R2) 
and an activation energy(E3) of 26.5 kcal/mol. An investigation of the relative stability of 
cyclic and open forms of the trimer by Karpfen 6 and coworkers has shown that the cyclic 
structure is more stable, while Liu et al. have investigated the stabilities of the trimer 
and tetramer relative to isolated fragments 7 . Most recently, Karpfen 8 has investigated 
equilibrium structures and concerted hydrogen exchange in (HF )3 and several other HF 
clusters, using split-valence plus polarization basis sets and the averaged coupled-pair 
functional (ACPF) method. These results are compared to ours below. 

In the present work we use large segmented Gaussian basis sets to locate the (HF )3 
stationary points at the SCF level, and then refine the energetics by performing extensive 
correlated calculations at the configuration-interaction and coupled- cluster levels, using 
large atomic natural orbital basis sets. Our computational methods are described in the 
next section, followed by a discussion of our results and conclusions. 


Computational methods 

The geometry of each stationary point was initially optimized at the SCF level using 
a triple- zeta plus polarization (TZP) basis of the form [5s 3p ld]/[3s 1 p] contracted from 
a (10s 6 p 2d/5s lp) primitive set 9 . The fluorine d set is a two-term fit to a Slater function 
with exponent 2.2 and the hydrogen p set has a Gaussian exponent of 1.0. Vibrational 
frequencies were evaluated at the SCF level using analytic second derivative methods . 
The TZP basis was extended to [5s 3p 2d] /[3s 2p] for use in MP2 calculations: an additional 
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d set was added (a two-term fit to a Slater function with exponent 0.7), while the original 
hydrogen p set was replaced by two functions with Gaussian exponents 1.4 and 0.35. 
Finally, in the largest of these MP2 calculations the [5s 3 p 2d]/[3s 2p] basis was augmented 
with an / set (exponent 1.9) on each of the fluorine atoms. All of these initial SCF and 
MP2 calculations utilized the program GRADSCF. 11 

For more elaborate studies of correlation effects, an atomic natural orbital (ANO) 
basis 12 of the form [4s 3p 2d If /3s 2 p Id] was used. This was contracted from a 
(13s 8p 6d 4//8s 6p 4 d) primitive set: the fluorine s and p exponents and the hydro- 
gen s exponents are taken from van Duijeneveldt 13 and the polarization functions are 
even-tempered expansions a(3 k , 0 < k < n. The ratio /? is 2.5 in all cases, with a(d) = 0.16 
and <*(/) = 0.49 on fluorine and a(p) = 0.1 and a(d) = 0.26 on hydrogen. The ANOs 
for fluorine were obtained from a single and double excitation Cl (CISD) calculation on 
the atomic ground state, while those for hydrogen were obtained from a calculation on the 

molecule H2. 12 

The first set of infinite-order correlated calculations was performed using the CISD 
method, including also Davidson’s correction for higher excitations 14 , denoted CISD+Q. 
The second method used was the coupled-pair functional (CPF) method of Ahlnchs and 
coworkers 15 which is nearly size-extensive. Finally, the coupled cluster method with single 
and double excitations (CCSD) was used — this is exactly size extensive 16 — together with 
the method denoted CCSD(T) in which a perturbational estimate of the effect of connected 
triple excitations is included. 17 In the correlated calculations SCF orbitals were used, and 
either 24 or 18 electrons were correlated: the former corresponds to neglecting fluorine 
Is electron correlation while in the latter correlation is neglected for both fluorine Is and 2s. 
No virtual orbitals were deleted in any calculations, and in the ANO basis calculations only 
the spherical harmonic components of the basis functions were used. The CISD and CPF 
calculations were performed using the MOLECULE-SWEDEN suite of programs 18 , while 
the coupled-cluster calculations were performed using the program VCCSD 19 . 

Results and Discussion 

Geometry 

Our optimized SCF bond distance in HF is 0.902 A, compared to an experimental 20 
value of 0.917 A. In the calculated equilibrium structure for the trimer, the H-F bond 
increases slightly, to 0.910 A, and the hydrogen- bonded H-F distance is 1.917 A. Our 
calculated results are similar to those of Gaw and co- workers 4 , who found r(H-F) to be 
0.912 A in the monomer and 0.923 A in the C 3h trimer. Their value for the H-F hydrogen 
bond distance is 1.860 A, somewhat shorter than our value, which was determined using 
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a larger basis set. Our calculated H-F bond distance in the symmetric D 3h trimer is 
1.136 A, in excellent agreement with the value of 1.139 A found by Gaw et al. As we noted 
above, Karpfen 8 , has included correlation at the ACPF level using a split-valence plus 
polarization basis set. His optimum bond lengths are 0.919 A in the monomer and 0.932 A 
in the trimer. Correlation thus appears to have a rather small effect on the geometries, 
although the effect is larger on the trimer than on the monomer. 


Vibrations: 

For HF we calculate a harmonic frequency of 4455 cm -1 with an infrared intensity 
of 172 km/mol using the TZP basis set at the SCF level. The experimental value is 
4138 cm -1 , while Gaw et al. calculated a value of 4440 cm -1 , and Karpfen 8 , a value of 
4182 cm -1 at the ACPF level. There is a strong dependence of the calculated vibrational 
frequency on the HF bond length: if the bond is elongated the computed “frequency” 
decreases sharply. The calculated trimer frequencies are given in Table I. 

The hydrogen- bonded systems (HF)„ have been extensively investigated experimen- 
tally, Some of these studies have included investigations of the higher species, although 
most have focused primarily on the dimer. Klemperer and coworkers 21 have shown that 
the trimer of HF is non-polar, consistent with a C 3h symmetry structure. Molecular beam 
predissociation experiments on the trimer have also been performed 22 ’ 23 , as have infrared 
vibrational spectroscopy experiments in neon and argon matrices 24 ’ 25 . The vibrational 
predissociation spectrum of the trimer indicates that the e' band lies at 3712 cm 1 . Only 
one of the bands observed in matrix isolation spectra is consistent with a Czh symmetry 
structure, the others are presumably due to open-chain forms that may be present in these 
experiments. This band is found at 3706 cm -1 in a neon matrix, and at 3702 cm -1 in an 
argon matrix. 

Our Cih structure is a minimum, based on the computed force constants. Of the 
twelve possible vibrational bands in this molecule there are three infrared active degenerate 
e' bands and one infrared active a" band. Two of these bands, one of a and the other 
of e 1 symmetry, correlate with the separated HF molecule fragment vibrations, while the 
remaining bands in the trimer correlate with translations and rotations in the separated 
fragments, that is to “intermolecular” modes in the trimer. The symmetric HF stretch is 
red shifted by 241 cm -1 while the degenerate HF stretch is red shifted by 175 cm -1 from 
the monomer. The red shift is consistent with the increase in the HF bond distance on 
going from monomer to trimer, although other factors also play a role, of course. Similar 
red shifts were computed by Karpfen. Our calculated frequencies are somewhat larger 
than those of Gaw and coworkers. The intensity of the degenerate e stretch is predicted 
to be very large both in our calculations and in the other studies. 
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The remaining harmonic frequencies are all below 1000 cm 1 and are associated with 
the intermodular librational motions. The symmetric and degenerate hydrogen bond 
bends are at 879 and 531 cm” 1 , respectively. The degenerate hydrogen bond bend is 
predicted to be very intense as is the out-of-plane motion of a symmetry. 

In order to compare our computed trimer frequencies with experiment, we scale them 
by a factor of 0.93, which is the ratio of the experimental harmonic frequency of HF 
monomer to our computed value. Our scaled calculations thus predict a harmonic fre- 
quency of 3919 cm -1 for the symmetric H-F stretch and 3980 cm 1 for the degenerate 
H-F stretch. The experimental 23 , gas phase result is 3712 cm -1 differing from our scaled 
result by almost 270 cm" 1 , and suggesting a remarkably large anharmonicity for this mode. 
The assignment of the remaining bands appears to be in much better agreement with the 
available experimental data. We assign the band at 590 cm 1 to the a out of plane mode, 
calculated value of 596 cm" 1 , while the band at 477 cm" 1 is assigned as the degenerate 
bend in the C 3h trimer, calculated value 478 cm" 1 . Both of these modes are predicted to 
be intense, while the lowest lying e' band is calculated to be very weak when compared to 
the other transitions. 

Force constant calculations confirm that the D^h symmetry structure is a transition 
state with the direction of negative curvature appropriate for the simultaneous exchange 
of hydrogens among the fluorines. The magnitude of the imaginary frequency is found 
to be 2224i cm" 1 , consistent with the dominant motion involving hydrogens. The totally 
symmetric stretch is about half of its value when compared to the Czh symmetry trimer. 
One of the degenerate stretches is much lower (1674 cm" 1 ) than the value found in the 
C^h symmetry trimer (4280 cm" 1 ), and one is much higher, 1506 cm versus 531 cm 

Energetics: 

Our discussion of the energetics of the concerted hydrogen exchange is based on the 
three reactions R1-R3 given in the Introduction: we repeat them here for convenience. 


3HF — *■ (HF) 3 [C 3h ] 

(Rl) 

(HF) 3 [C 3h ]^(HF) 3 [D 3h ] 

(R2) 

3HF — ► (HF) 3 [D 3k ]. 

(R 3) 


Although ultimately the enthalpy changes are important chemically, we will first focus 
only on the electronic energy differences. Our results are summarized in Table II, where 
we present values for each of the three reactions given above with SCF, perturbation theory, 
Cl, as well as the coupled-cluster methods. Our most accurate results suggest that the 
stable hydrogen bonded trimer is bound by slightly more than 14 kcal/mol relative to the 
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separated fragments, the barrier to reaction (H3) is predicted to be only 3.6 kcal/mol, 
and the energy separation between the transition state, and the stable trimer (-R2) is 
18 kcal/mol. This is in marked contrast to the results of most of the previous work cited 
in the introduction. There is also a substantial variation in the results depending on the 
correlation treatment employed, as we now discuss. 

At the CISD+Q level of correlation treatment, the energy difference between three 
separated HF molecules and the transition state is found to be small: about 7 kcal/mol. A 
comparison of this reaction with the covalent exchange of hydrogens in the H 6 system indi- 
cates that introduction of ionic character into the molecular system by replacement of three 
hydrogens with three fluorine atoms reduces the effective barrier from some 70 kcal/mol 
to about 7 kcal/mol. An approximate quantitative measure of these changes can also be 
seen by an examination of the Mulliken charges calculated for the trimer and the transi- 
tion state. At equilibrium the charge on hydrogen is found to be 0.3e _ , whereas at the 
transition state this increases to 0.4e~. 

As expected, the formation of the C 3h trimer from three monomers is found to be 
exothermic. The SCF value is found to be 12.4 kcal/mol while each of the correlated 
treatments predict that the complex is slightly more than 14 kcal/mol more stable than 
the reactants. The effect of correlation is found to be about 1.7 kcal/mol and roughly 
independent of the number of electrons correlated or the level of treatment. Liu and 
coworkers ‘ have found that correlation contributed about 1 kcal/mol to this energy, some- 
what less than we do. Their value was obtained with a TZP basis, while correlation effects 
were included at the ACCD level^ . Karpfen obtains a somewhat larger binding energy 
(15.4 kcal/mol) at the ACPF level, with a rather small basis. We note, however, that 
consideration of basis set superposition error changes these observations, as we discuss in 
detail below. 

The effect of correlation is most dramatically seen on the energy difference between 
the stable complex and the transition state. At the SCF level we find this energy differ- 
ence to be over 35 kcal/mol, while the simplest correlation treatment (MP2) reduces this 
energy difference to 16 kcal/mol. More elaborate single and double excitation correlation 
treatments result in a difference of about 20 kcal/mol. The results are largely insensi- 
tive to whether 18 or 24 electrons are correlated, so it would not be unreasonable here to 
regard the fluorine 2s electrons to be part of an atomic core. On the other hand, there 
are significant differences between the CISD results and those of the various size-extensive 
or approximately size-extensive treatments; with such a large number of electrons corre- 
lated this would be expected. Thus the simplest Cl treatment, CISD, predicts an energy 
separation of 23.9 kcal/mol, whereas inclusion of the +Q correction reduces this value to 
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21.1 kcal/mol. CPF reduces this separation by another 1.1 kcal/mol, but this appears to 
be an overestimate when compared to the CCSD result of 20.9 kcal/mol. Finally, inclu- 
sion of the (T) correction for the effects of connected triple excitations reduces the energy 
separation by almost 3 kcal/mol. This large triples contribution to the barrier height for 
hydrogen exchange deserves further comment. 

On initial examination, such a large effect of triples is unexpected, and might be 
assumed to indicate that nondynamical correlation (near-degeneracy effects) might be im- 
portant at the transition state. Recently, Lee et al. 27 have proposed the use of a diagnostic 
for the effects of nondynamical correlation, based on the magnitude of the norm of the 
amplitudes, denoted T x . This T x diagnostic has been shown to be a good indicator of 
the applicability of a single-reference-based coupled- cluster method. It was suggested 
that a value larger than 0.02 for 7) indicates nondynamical correlation effects are large 
enough to cast doubt on the reliability of single-reference treatments limited to single and 
double excitations. For the D 3k geometry of (HF) 3 T x is 0.01, which would indicate that 
nondynamical correlation effects should not be a problem. The large triples contribution 
in fact arises from a different (and in some respects a simpler) cause. At the transition 
state we have three chemical bonds which are significantly elongated relative to their equi- 
librium value. We can examine the effect of triple excitations on this bond-lengthening by 
computing the triples contribution for a single HF moiety, at equilibrium bond length and 
at the bond length obtained for the transition state. The triples contribution to the energy 
difference between these two geometries was found to be a little less than 0 . 1 kcal/mol. 
Thus of the 2.9 kcal/mol triples correction to the barrier height, almost 2.1 kcal/mol origi- 
nates in the simultaneous stretching of three HF bonds. Hence only 0.8 kcal/mol should be 
ascribed to a true “many-body” effect. Nevertheless, it is clear that in a system where even 
more bonds participate in the rearrangement, such as hydrogen exchange among higher 
oligomers of HF, such cumulative effects would be further magnified and would certainly 
require an appropriate correlation treatment. 

We can now examine several sources of error in these calculations. The first is the 
use of SCF geometries in the correlation treatment. In order to address this question 
we have reoptimized the structure of the D^h transition state at the MP2 level using the 
TZ2P basis augmented by a set of /-type functions on each of the fluorine atoms. The 
bond lengths increase slightly (less than 0.01 A), while the energy is lowered by less than 
0.8 kcal/mol relative to the MP2 result at the SCF geometry. Although we have performed 
only the easier task of reoptimizing the geometry of the D^h symmetry structure with its 
two degrees of freedom, we expect the effects of electron correlation to be substantially 
less for the equilibrium structure than for the transition state, as discussed below. We are 
therefore confident that our choice of SCF geometries will not be a source of significant 
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error. 

A quantitative assessment of our energetics must include an estimate of the basis set 
superposition error (BSSE). Ideally, this should be done for both the stable complex and the 
transition state. However, estimating the BSSE for polyatomic systems like (HF) 3 presents 
some conceptual difficulties. We have chosen to proceed by computing a counterpoise 
correction for the stable complex, using the SCF-optimized geometry and with both a 
single HF molecule with two HF ghost basis sets present and two HF molecules with 
a single ghost basis present as the fragments. The ANO basis was used in all BSSE 
investigations. At the SCF level the BSSE is less than 0.1 kcal/mol. However, at the 
CCSD level the BSSE is slightly greater than 0.5 kcal/mol for each HF moiety. We would 
thus assign a BSSE of 1.5 kcal/mol for the hydrogen-bonded trimer: this is likely to be 
an overestimate rather than an underestimate. Nevertheless, these results suggest that 
the effects of correlation on the binding of the trimer are grossly exaggerated, since most 
of the energy lowering is due to BSSE. This is consistent with the results of Karpfen, 
who (assuming an SCF result similar to ours) would obtain an even greater correlation 
contribution to binding: the split-valence plus polarization basis of Ref. 8 is likely to have 
a larger BSSE than our ANO basis. On the other hand, the results of Liu and co-workers 
are less consistent with this picture. Finally, the structure of the transition state is not 
qualitatively different from that of the stable complex, so we can reasonably assign the 
same (upper bound) value of 1.5 kcal/mol to the BSSE in the transition state. 

The incompleteness of the one-particle space is also a source of error in these calcula- 
tions. Comparison of the MP2 results, at least, in the TZ2P and ANO basis sets, suggests 
that the effect of basis set extension on the binding energy of the trimer is relatively small. 
Nevertheless, it is now well understood that the dominant contribution to hydrogen bond- 
ing is electrostatic, and so the basis set used should not only correctly describe the bonding 
in the monomer, but also provide a reasonable description of the multipole moments and 
polarizability. ANO basis sets are not always capable of meeting these requirements with- 
out further augmentation 28 , which would make the present calculations too expensive. It 
is thus possible that the effect of basis set incompleteness on the binding energy is some- 
what underestimated by our calculations. However, it seems unlikely that the remaining 
effects would contribute more than 1 kcal/mol. Since the bonding in the transition state 
is similar to that in the stable complex, we expect that the basis set effect on the barrier 
height would, if anything, be even smaller. 

Enthalpy Changes 

Our calculated A E values refer to the bottom of the potential energy well and need 
to be corrected by the difference in zero point vibrational energy, AZPE. Our correction 
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is done at the harmonic level using the experimental value for the HF frequency, while 
the trimer values are derived from our calculated frequencies scaled by 0.93 as indicated 
above. The stable trimer contains 4.5 kcal/mol more vibrational energy than the separated 
reactants which reduces its stability to 9.9 kcal/mol. Correcting this value for BSSE will 
further reduce the stability of the trimer: if the full counterpoise correction is applied the 
result would be as low as 8 kcal/mol. In view of our discussion of basis set incompleteness, 
it seems reasonable to assert that the trimer is bound by 9il kcal/mol. 

The AZPE correction between the transition state and stable trimer is -3.5 kcal/mol, 
while it is 1.1 kcal/mol between the transition state and the reactants. Given these values 
our calculated barrier at 0 K is 4.7 kcal/mol, and the energy separation between the stable 
trimer and the transition state is 14.5 kcal/mol. Each of these values may then be further 
corrected for BSSE, which would raise the computed barrier to 6.2 kcal/mol. 

These results require further correction for non-zero temperatures. For reaction R2 
there are no corrections required to convert A E to A H if we neglect the dependence of 
the vibrational energy on temperature. For R1 and R3 the values must be corrected for 
differences in the translational and rotational energies and for the difference between A E 
and AH (AH = A E + AnRT with An = -2). The correction term is -6.5 RT and 
at 300 K this corresponds to a value of — 3.9 kcal/mol for the correction, or a barrier of 
2.3 kcal/mol, if the BSSE value for the barrier is used. 

Conclusions 

We have investigated the exchange of hydrogen atoms in a hydrogen-bonded sys- 
tem and have shown that such an exchange can proceed with an activation energy of 
~4 kcal/mol. In addressing this problem we have also demonstrated that the presence of 
ionic character in the bonding lowers the barrier relative to that found in a purely cova- 
lent system, such as H 6 . The effects of electron correlation on the energetics once again 
demonstrate the magnitude of size-extensivity errors when many electrons are correlated. 
In addition, in (HF) 3 the effects of triple excitations significantly alter the barrier height 
from that predicted at the CCSD level, and we can expect this effect to become more im- 
portant with increasing size of the system. Finally, since the barrier to concerted hydrogen 
exchange is small, we may anticipate a significant contribution from tunneling to the rate 
of the exchange reaction. 
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Table I 

Calculated harmonic vibrational spectra 


Mode uj (cm -1 ) I(km/mol) Description 



C^h Structure: 



a 7 

4214 

0 

symmetric HF stretch 

a 7 

879 

0 

symmetric bend 

a 7 

186 

0 

symmetric stretch -libration 

a" 

662 

507 

out of plane torsion 

e 7 

4280 

571 

degenerate HF stretch 

e' 

531 

382 

degenerate angular def. 

e' 

157 

17 

degenerate stretch - libration 

e" 

475 

0 

degenerate out of plane torsion 


Dzh Structure: 



a '\ 

2212 

0 

symmetric breathing stretch 

a\ 

788 

0 

symmetric HF str. and bend 

a 2 

2194i 

- 

negative curvature stretch 

a'± 

1449 

466 

out of plane torsion 

e< 

1674 

21 

degenerate stretch plus bend 

e' 

1506 

4154 

degenerate stretch + bend 

e' 

604 

0.3 

degenerate bend + stretch 

e" 

1135 

0 

degenerate out of plane torsion 
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Table II 

Comparison of relative energies" 


Calculation: 

AF(R1) 

A£(R2) 

A£(R3) 

SCF, Cl, and CPF Results: 



ANO SCF 

-12.4 

35.8 

23.4 

ANO CISD 18 el 

-14.2 

23.8 

9.6 

ANO CISD+Q 18 el 

-14.1 

21.3 

7.2 

ANO CPF 18 el 

-14.1 

20.8 

6.7 

ANO CISD 24 el 

-14.3 

23.9 

9.6 

ANO CISD+Q 24 el 

-14.2 

21.1 

6.9 

ANO CPF 24 el 

-14.1 

20.0 

5.9 

SCF, MP2, and CC Results: 



TZ2P/SCF 


35.1 


TZ2P/MP2 


16.7 


ANO SCF 

-12.4 

35.8 

23.4 

ANO/MP2 

-14.3 

16.1 

1.8 

ANO CCSD 24 el 

-13.9 

20.9 

7.0 

ANO CCSD(T) 24 el 

-14.4 

18.0 

3.6 

a) Relative energies are given 

in kcal/mol. 
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